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The  objective  of  this  project  is  to  combine  a  number  of  recent 


advances  in  finite  element  theory  and  computer  technology  for  analyzing 
cavities  and  structures  in  rock.  This  computer  program  applies  to  general 
three-dimensional  structures,  considers  nonlinear  material  properties  includ¬ 
ing  homogeneous  deformations  and  inhomogeneous  deformations  due  to  joints, 
anisotropic  and  time-dependent  material  propertoes,  gravity  loading,  and 


sequence  of  construction  or  excavation. 


During  the  first  half  of  this  contract,  work  has  been  aimed  at  pro¬ 
ducing  a  user-oriented  computer  program.  The  work  of  writing  the  program  was 
divided  into  three  areas: 


a.  Input 

b.  Execution  and  output 

c.  Materia]  properties 


The  Input  Section  automatically  generates  the  continuum  part  of  the  finite 
element  mesh,  including  joint  elements,  allows  the  user  to  add  other  elements 
(beam,  shell,  truss)  to  the  mesh,  plots  the  result,  reduces  the  bandwidth  and 
reads  loads,  material  properties,  and  other  quantities  necessary  to  the  cal¬ 
culation.  The  Execution  Section  forms  the  global  stiffness  matrix  and  solves 
equations  of  equilibrium  for  displacements  by  an  implicit  method.  The  material 
properties  are  represented  by  subroutines  within  the  Execution  Section,  which 
are  written  in  a  modular  form  so  that  if  the  general  equations  of  nonlinear 


elasticity,  viscoelasticity,  viscoplasticity,  or  plasticity  do  not  suit  a 


particular  problem  they  may  be  easily  modified.  This  work  is  now  complete 
except  for  linkage  among  the  various  sections  to  be  checked  out  and  for  the 
efficiency  of  some  operations  to  be  improved. 
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One  of  the  guidelines  for  this  project  was  that  no  major  new  research 
was  to  be  done.  Accordingly,  the  program  uses  existing  finite  elements,  a 
proven  form  of  the  equation  of  equilibrium,  existing  material  property  descrip¬ 
tions,  and  an  existing  bandwidth  reducer.  However,  a  small  amount  of  new  work 
was  done.  A  new  joint  element  was  developed  and  an  existing  concept  for 
automatic  mesh  generation  was  greatly  extended.  Also,  a  form  of  Choleski 
decomposition  was  modified  for  efficient  use  of  multibuffering,  resulting  in 
substantial  improvement  in  efficiency  of  solving  equations  of  equilibrium  using 
peripheral  storage. 


During  the  first  six  months  of  this  study,  some  new  work  was  reported 
by  other  ARPA/Bureau  of  Mines  contractors  which  has  been  incorporated  in  the 
program.  Among  these  are  some  creep  data  obtained  by  W.  A.  Wawersik  of  the 
University  of  Utah  and  s trength/deformabl 1 1 ty  data  for  faults  by  R.  E.  Goodman, 

F.  E.  Heuze,  and  Y.  Ohnishl  of  the  University  of  California,  Berkeley. 

The  technical  work  reported  below  is  divided  Into  four  main  parts. 
Section  3  describes  the  equations  of  equilibrium  and  the  types  of  finite  elements 
available  in  the  program,  including  the  new  joint  element.  Section  k  describes 
the  properties  of  rock  including  anisotropy,  rate  effects,  joints,  and  homogen¬ 
eous  inelasticity.  Section  5  describes  the  computer  program.  Most  attention  is 
given  to  the  processing  of  input  data  and  to  options  available  to  the  user  for 
mesh  generation,  sequential  excavation,  or  construction,  automatic  bandwidth 
reduction  and  plotting  of  the  mesh.  Some  logical  flow  diagrams  are  given. 

Section  6  describes  some  example  problems  which  have  been  solved  during  the 
checkout  of  the  program. 

During  the  second  half  of  this  project  work  will  concentrate  on 
applications  of  the  program  to  field  situations  In  which  measurements  are 
available  for  comparison  with  results  of  analysis. 
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SECTION  1 
INTRODUCTION 

The  purpose  of  this  contract  is  to  combine  a  number  of  recent 
a  Vances  finite  element  theory  and  computer  technology  into  a  computer 
program  or  analyalng  structures  and  cavities  in  rock.  This  computer  program 
appl les  to  general  three-dimensional  structures,  considers  nonlinear  mate rial 

,n:'ud'n9-horo9eneous  deforrau°"  a-d  -.*~.ti.n  due 

.  an  sotropic  and  time-dependent  material  properties,  gravity  loading 

and  sequence  of  construction  or  excavation  finr  m  9* 

DMrf.  ,  ,  ,  excavation.  u,nce  the  program  is  intended  for 

cultles  ^  ^  and  deSi9n’  9re3t  eff°rt  1,35  6ee"  ™ada  “  foresee  dlffi- 

bv  th  T9  ''  F°r  eXa,"p,e'  much  tedious  work,  which  formerly  was  done 

width*  T’  bee"  e""inated  by  s°Ph'”'oated  mesh  generators  and  a  band- 
re  ucer.  Also,  since  many  prospective  users  may  have  access  to  small 

000  Z5'  “mP:tefS  Sti"  Wi5b  “  Problems 

(4000  6000  eguations),  the  program  uses  up-to-date  mu, t ibuf fer ing  technigues 

^accessing  peripheral  storage  units,  thus  dramatically  reducing  compel  r  run 
use"  M,r  Til0"  Pr06'emS-  Fina"V'  ^  aUemPt  '*  “  '«.«•■»  the 

fo  1  *  mate.r,,a'  Pr°Perty  deSCriPtio"  and  by  making  the  program  eff  iclent 
for  and  compatible  with  a  wide  variety  of  computers. 

The  purpose  of  this  report  is  to  discuss  what  the  computer  program 
ow  oes  and  what  it  will  do  when  the  present  contract  is  complete.  Most  of 
-script, on  ,s  given  from  the  standpoint  of  a  prospective  user.  Mesh 
ganeretion  the  types  of  elements  available,  the  types  of  loading  and  construction 

ZtZ  a'0"'  and  r  . . .  aPd  J°'ntS  ”hlCh  a™  ayailabla  are 

analv  ,  T  a‘  Pr°9ra"  "  W°rki"9  pr°parly  is  ■  “v  -"paring 

element  sT  !  ”  ^  °ne"'  ^ ^  three‘d  '“a"5  Iona!  problems  with  the  finite 
n  solutions  to  the  same  problems.  The  theory  underlying  the  present  finite 

ement  formulation  is  also  described.  Also,  the  structure  of  the  code  is 
indicated  by  logic  diagrams. 
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The  work  which  will  be  done  in  the  second  half  of  the  contract  period 
involves  application  of  the  code  to  actual  mining  situations.  During  this 
phase,  general  features  of  the  program  such  as  the  joint  element  and  the 
material  property  representations  will  be  refined.  Parameters  will  be 
selected  for  specific  field  conditions  and  an  attempt  will  be  made  to  match 
field  data  with  calculations.  The  sites  have  not  yet  been  selected. 

To  summarize,  this  report  describes  a  finite  element  computer  program 
which  incorporates  a  number  of  good  features  which  are  not  all  available  in  any 
other  program.  By  means  of  comparison  with  analytic  solutions,  evidence  is  given 
that  the  program  is  consistent  with  the  assumptions  under  which  it  was  formu¬ 
lated  and  written,  and  that  It  contains  no  programming  errors.  The  task  of 
demonstrating  that  it  accurately  represents  field  conditions  is  left  to  the 
second  phase  of  the  contract. 
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PHYSICAL  ASPECTS  OF  ROCK  AND  SUPPORT  SYSTEMS 
WHICH  THIS  PROGRAM  REPRESENTS 

This  section  describes  the  scope  and  specific  goals  of  the  prWnt 
computer  program.  These  are  representation  of  gravity  and  live  loading  v 
properties  of  homogeneous  rock  including  anisotropic,  inelastic  and  Viscous 

e  fects,  properties  of  jointed  rock,  excavation  and  construction,  and'  geometry 
of  cavities  and  support  systems. 

To  illustrate  the  capability  of  the  program  three  typical  problems 
are  described  below.  These  problems  have  not  actually  been  solved  with  .the 
present  program,  but  they  could  be  solved  at  any  time.  The  first  Is  Ulus- 
trated  in  Figure  2-1.  A  section  of  tunnel  is  to  be  excavated  in  a  region 
containing  a  major  Joint.  The  properties  of  the  Joint  are  assumed  to  be  known. 
The  rock  adjacent  to  the  joint  is  assumed  to  be  homogeneous  and  to  have  viscous 
properties  which  can  be  represented  by  a  visco-plastic  model.  In  Step  I,  the 
tunnel  has  not  yet  been  excavated.  Stress  in  the  rock  Is  computed  by  applying 
static  overburden  to  the  edges  of  the  finite  element  mesh.  Then  the  tunnel  is 
excavated  by  removing  appropriate  elements.  At  each  stage  of  the  excavation 
the  tunnel  roof  Is  propped  by  truss  and  beam  elements.  Eventually,  the 
tunnel  Is  fully  open  and  the  final  supports  are  Installed.  Each  stage  Is 
associated  with  an  elapsed  time,  during  which  the  rock  flows  In  a  vlsco- 
plastlc  manner.  At  each  intermediate  stage  and  at  a  stage  the  user  defines 
to  be  final,  the  stresses  in  the  rock  and  in  the  support  elements  are  printed. 

The  second  problem  is  Illustrated  In  Figure  2-2.  A  bank  Is  to  be 
excavated  In  a  rock  such  as  shale  having  nonlinear,  anisotropic  stress/straln 
properties  and  an  anisotropic  fracture  criterion.  In  Step  I,  the  In  situ  states 
o  stress  are  computed  by  applying  gravitational  forces  In  a  step-by-step 
fashion  throughout  the  grid.  In  subsequent  steps,  elements  are  reneved  In  any 
sequence  the  user  desires.  Between  excavation  steps  the  remaining  rock  will 

be  checked  for  fracture  which  would  correspond  to  spall  and  sliding  In  an 
actual  field  situation. 
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BEDDING  PLANE 


STEP  2.  EXCAVATION 


STEP  3.  EXCAVATION 


FIGURE  2-2.  EXCAVATION  OF  BANK 
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The  third  example  Is  an  extension  to  three  dimensions  of  the  first 

example.  The  final  stage  In  the  calculation,  at  which  the  section  of  tunnel 

under  consideration  is  fully  excavated,  is  illustrated  In  Figure  2-3.  Notice 

that  the  plane  of  the  joint  is  not  parallel  to  the  axis  of  the  tunnel.  This 

problem  is  one  of  the  largest  and  most  difficult  which  the  code  is  Intended  to 
handle. 


Several  Important  questions  remain  unanswered.  First,  how  much  do 
such  analyses  cost  In  terms  of  computing  time  and  engineering  preparation  and 
Interpretation?  Second,  are  sufficient  field  data  available  for  nonlinear  and 
anisotropic  properties  of  the  rock  and  for  the  properties  of  the  joints  to  be 
accurately  modeled?  Part  of  the  answer  to  the  first  question  can  be  found  In 
Section  6  of  the  present  report  In  which  the  computer  times  for  solving 
simple  checkout  problems  are  given.  More  information  on  costs  will  be  forth¬ 
coming  at  the  end  of  this  project.  It  appears  at  present  that  sufficient 
field  data  are  not  available  for  accurate  modeling.  Hopefully,  the  develop¬ 
ment  of  analytic  tools  such  as  the  present  program  will  stimulate  the 


necessary  measurements. 
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FIGURE  2-3.  THREE-DIMENSIONAL  TUNNEL  ANALYSIS  WITH  JOINTS 


SECTION  3 

APPLICATION  OF  FINITE  ELEMENT  THEORY 
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This  section  discusses  the  present  formulation  of  equations  of 
equilibrium.  The  provisions  to  extend  this  formulation  to  1/irge  deformations  Is 
also  described.  Then  the  typss  of  elements,  Including  truss,  beam,  plane  strain 
and  axlsymmetrlc,  three-dimensional  and  thick  shell  elements  are  described,  in 
addition,  a  new  element  for  representing  slip  and  debondtng  along  planar  joints 
is  described. 

3*1  SOLUTION  OF  NONLINEAR  EQUATIONS  OF  EQUILIBRIUM 

The  matrix  equation  of  equilibrium  for  a  structural  system  with 
material  nonlinearity  Is: 

K  (u)  u  -  P  (3_ | ) 

where  the  Instantaneous  stiffness  matrix  (K)  Is  a  nonlinear  function  of  the 
displacement  vector  (u) .  P  Is  the  vector  containing  external  loads.  There 
exist  numerous  methods  of  solving  the  above  system  of  nonlinear  equations.  In 
general,  these  methods  can  be  divided  into  two  classes:  Iterative  methods  and 
Incremental  methods. 

Iterative  methods  apply  the  total  load  Initially  and  approach  the 
solution  by  modifying  the  stiffness  matrix  and/or  modifying  the  load  vector. 
Modification  of  the  stiffness  matrix  In  general  accelerates  the  convergence 
but  is  computationally  costly.  The  optimum  may  be  achieved  by  occasional 
stiffness  reformulation; 
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In  incremental  methods  the  loads  are  applied  in  several  steps  and  an 
incremental  form  of  the  Equation  3-1  is  solved. 


where 


Au  . , 
—  n+1 


-n+l 


u 

-n 


P 

-n 


It  is  important  to  note  that  the  stiffness  matrix  can  only  be  formed  based 

on  the  displacement  vector  from  the  previous  step  u  which  creates  some 

-n 

step-wise  error,  in  tills  simple  incremental  technique  the  step-wise  errors 
can  accumulate  and  lead  to  considerable  total  error.  To  prevent  this  accumula¬ 
tion  of  the  step-wise  error,  a  modified  form  of  the  load  vector  Is  used. 


K  Au  , 
-n  — n+l 


(3-3) 


where 

En+j  "  T°tal  load  vector  at  the  end  of  the  (n+l)th  step 

-n  "  Vector  of  the  internal  resisting  forces  at  the  end  of  the 
th 

n  step 

By  using  this  method  of  load  vector  correction  the  equilibrium  is  satisfied  at 
the  beginning  of  each  Incremental  step  and  thereby  the  accumulation  of  the 
step-wise  error  Is  prevented.  Satisfaction  of  equilibrium  Is  assured  in  spite 
of  errors  or  approximations  In  the  stiffness  matrix  and,  therefore,  the  reformu 
lation  of  the  stiffness  matrix  is  not  required  at  every  step.  However,  the 
error  in  each  step  is  directly  dependent  on  the  approximation  of  Instantaneous 
stiffness  matrix. 
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An  alternative  method  is  to  apply  the  total  force  from  the  beginning, 
In  which  case  the  load  in  Equation  3-3  will  be 

Pn  "  P  total  (n“‘ . N> 

It  should  be  noted  that  the  application  of  the  total  loads  makes  this  method 
equivalent  to  an  Iterative  scheme  with  load  vector  correction  which  was 
previously  discussed.  However,  the  loads  In  general  have  a  specified  history 
dictated  by  the  sequence  of  application,  sequence  of  construction  and  excava¬ 
tion,  and  the  time  phenomenon  associated  with  the  viscous  material  properties. 

In  most  practical  problems,  the  specified  history  of  loading  is  a  series  of 
step  functions.  This  is  true  in  case  of  construction  and  excavation  which  can 
be  considered  as  a  discontinuity  In  force-displacement  relation  and  an  abrupt 
change  in  the  instantaneous  stiffness  matrix. 

An  efficient  scheme  is  to  apply  the  total  of  the  step-wise  loading 
at  each  stage  and  then  carry  out  several  Iterations  with  occassional  stiffness 
reformulation  to  accelerate  the  convergence.  This  scheme  is  summarized  in  the 
following  steps. 

For  each  step: 

a.  Compute  un  -  un_j  +  Aun  (for  first  step;  uq  -  0) 

b.  Compute  the  strains  (e  )  or  strain  Increments  (&e  )  using  the 

n  -n 

derivatives  of  the  shape  functions  for  each  element  which  have 
been  Initially  computed  and  stored 

c*  For  time- independent  materials  compute  the  stress  (0  )  and 

the  Instantaneous  stress-strain  relations  (C  )  (see  section  on 

-n 

material  properties) 

2.  For  visco-elastic  elements 

1.  Compute  stresses  2n  -  C  (en  -  cjj,,)  where  C  is 

the  elastic  stress-strain  matrix,  e  |s  the  total  strain 

C  “n 

an<*  £n-j  is  the  total  creep  strain 
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li.  Using  the  stresses,  compute  e  which  is  the  total 

-n 

creep  strain  at  the  end  of  the  time  step  (see  section  on 

material  properties) 

111.  Compute  effective  stress 

o_  -  C  (e  -  cC) 

-n  -  '-n  -n 

d*  Compute  the  internal  resisting  forces  from  the  stresses  (effec¬ 
tive  stresses  for  viscoelastic  elements),  if  it  is  stiffness 
update  cycle,  compute  stiffness  matrix. 

e.  Solve  Equation  3_3  to  obtain  Aun+j .  Compute  y  «  ||  Au  ||„ 

f*  If  a  specified  number  of  iteration  has  been  reached  or  if 

y  <  e  (e  is  a  specified  quantity),  go  to  Step  g\  otherwise, 
go  to  Step  a  and  repeat  the  iteration. 

9*  Apply  the  next  loading  step  and  go  to  Step  a. 


3-2  EQUATIONS  OF  EQUILIBRIUM  FOR  LARGE  DEFORMATIONS 


The  method  of  large  deformation  finite  element  analysis  to  be  used  In 
the  present  computer  program  was  initially  introduced  by  Sharifi  and  Yates, 
Reference  3"1 . 


The  matrix  equations  of  equilibrium  (or  motion)  are  derived  from  an 
incremental  virtual  work  expression  and  the  original  configuration  of  the  finite 
element  system  is  taken  as  the  reference  configuration.  This  choice  of  the 
reference  state  eliminates  the  need  for  updating  of  the  coordinates  of  the 
nodal  points  which  is  computationally  a  costly  operation. 
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■  ion  Is 


nta,  virtual  work  express. 

The  Incremental 

c  s..  M*sfa  '  f,  S,J  S‘"iJ<IVo 

nv^-.'V/O1  U 

Jao  «-« 

+  f  as,,  sae'j  dV0 
•Ai  J 


where 


of  displacement  vector 

-  Component  displacement  vector 

*  .  incremental  component  of 

k  dant  of  traction  vector 

*  a,  component  of  traction  vector 

»  Incremental  comp 


&u 

t 


At. 


S 


U 


e.  Plola  stress 

„  Component  o 


tensor 


tensor 


ta,  component  of  Plo'e  *tre5S 

-  Increments  n  tenaor 

«nent  of  incremental 

.  Linear  componen  stra|h  tensor 

^  ij  _*■  of  incremental 

Ln  .  Monllnear  component  o 


*Slj 


of  the  or 


j  the  coordinates 

U  .  are  referred  to  tire  area  and 

Th.  stresses  and  tract, on  are 

,9,na|  conf l9urat'on:  (3-5a) 

+  uv  iauit  j  +  ulc* 
au,  ,  +  AUi.i  k-’  k,J 

“lJ  ,J  (3-5b) 


an 


U 


displacements  «!*■»  «‘h 

The  total  and  values  of  the  displacement 

shape  functions 


be 
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Ui  =  HU, 


(3-6a) 


Au j  =  H  AU. 


(3-6M 


where  H  is  the  Vector  of  the  shape  functions  and  U.  and  AU..  are  the 
vectors  of  the  total  and  incremental  displacements  of  the  element  nodes. 

Without  loss  of  generality  the  remaining  part  of  this  section  will 
be  devoted  to  the  derivation  of  the  appropriate  matrices  for  the  two- 
dimensional  quadrilateral  element. 

Substituting  the  Equation  3-6  into  Equation  3~5  will  result  in 
the  following  expressions  for  the  strain  increments  in  terms  of  the  nodal 
point  displacements: 

iExx  '  -’x  4HX  +  <S-x  V  <B-x  ^x’  ♦  <8'X  -y>  <B.x  ^ 


=  — ’y  4yy  +  <d.v  yx)  <d,y  4Ux)  +  (H,v  Uv)  (H,v  AU^) 


Ae  -  H,  AU  +  H,  AU  +  (H  U  )  (H,  AU  )  +  (H.  U  ) 

xv  -  x  -y  -  y  -x  -y  -x  -’v  » u 


y  -x' 


(3-7) 


(H,x  AUy)  *  (H,x  Ux)  (H,v  AUx)  +  (H>u  AUJ 


x  -x  -’y  -x 


y  -y 


The  above  equations  can  be  written  as  matrix  form 


Ae  =  B  (I  +  E)  AU  =  §  AU 


(3-8) 


where  I  is  the  identity  matrix  and 
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(3-9) 


H»x  -/ 

H,y  Uy  H,y 


(3-10) 


B  is  the  usual  strain-displacement  matrix  for  infinitesimal  deforma 
tion  and  E  is  the  large  deformation  contribution. 


The  linear  and  geometric  stiffness  matrices  and  the  load  correction 

vector  are 


be 


B  dv 

m 


B  dv 

m 


(3-11) 


(3-12) 


(3-13) 


where  v  is  the  area  of  each  element  in  original  configuration,  and  C  is 
m 

the  instantaneous  stress-strain  relation 


AS  =  C  A  e 

Finally,  the  matrix  equat'on  of  equilibrium  is 
(be  +  bg)  AU  =  R  -  F 


(3-H) 


(3-15) 
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It  is  important  to  note  that  for  the  computation  of  the  above  matrices 
only  the  derivatives  of  the  shape  functions  H,x  and  H ,y  at  original  geometry 
of  each  element  is  required.  Therefore,  these  derivatives  at  integration 
points,  can  be  computed  in  the  first  part  of  the  program. 

3.3  STRUCTURAL  FINITE  ELEMENTS 

A  description  of  the  structural  elements  incorporated  in  this  computer 
program  are  given  here.  The  beam  and  thick  shell  elements  have  linear  elastic 
properties.  All  other  elements  are  capable  of  representing  nonlinear  properties. 

3.3.1  THREE-DIMENSIONAL  TRUSS  ELEMENTS 

The  truss  element  is  the  conventional  space  truss  member  which  can 
resist  compression  or  tension  along  its  axis.  It  can  also  be  used  to  model 
bolts.  The  truss  member  is  subject  to  three  translations  at  each  end  of  the 
member  as  shown  in  Figure  3-la.  The  member  stiffness  matrix  is  of  order 
6x6.  The  material  and  geometrical  properties  are  defined  by  the  tangent 
Young's  modulus,  and  the  cross-sectional  area  of  the  element. 

3.3.2  THREE-DIMENSIONAL  BEAM  ELEMENTS 

The  three-dimensional  beam  element  is  subject  to  three  translations 
and  three  rotations  at  each  end  of  the  member.  The  generalized  forces  and  the 
generalized  displacements  associated  with  the  s ix-degrees-of-freedom  (DOF)  at 
each  end  are  shown  in  Figure  3-lb. 

The  geometrical  properties  of  the  beam  element  are  specified  by  an 
axial  and  two  shear  areas  and  three  principal  moments  of  inertia,  two  associated 
with  bending  and  one  with  torsion.  Young's  modulus  and  Poisson's  ratio  are 
required  to  define  the  material  properties  of  the  beam  element. 

The  element  stiffness  matrix  is  of  order  12  x  12  and  is  obtained 
from  the  classical  beam  theory  including  the  effects  of  the  shear  deformations. 
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(a)  THE  THREE-DIMENSIONAL  TRUSS  ELEMENT 


(b)  THE  THREE-DIMENSIONAL  BEAM  ELEMENT. 
FIGURE  3-1.  TRUSS  AND  BEAM  ELEMENT 


16 


R-721 5-2299 


A 


A  provision  for  the  member  end  boundary  conditions  accounts  for 
hinges  and  other  releases. 

3.3.3  TWO-DIMENSIONAL  PLANE  STRAIN  AND  AXISYMMETRIC  ELEMENTS 

Quadrilateral  Isoparametric  elements  will  be  used  in  the  computer 
program.  For  a  qenoral  quadrilateral  element,  as  shown  in  Figure  3-2,  the 
local  and  global  coordinate  systems  are  related  by 


X  3  Z  h.  x. 


y  =  A  K  y* 


where  the  Interpolation  functions  are  given  by 


(3-16) 


hl  3  1/4  (1-s)  (1-t) 
h2  *  1/4  0+s)  (1-t) 
h3  55  1/4  (1+s)  (1+t) 
h4  3  1/4  (1-s)  (7n) 

The  same  Interpolation  functions  are  used  in  the  displacement 
approximation. 

f 

UX  *  Z  h1  Ux1  +  h5  a|  +  h6  a2 

Uy  (S.t)  =  Z  h1  Uy1  ♦  h5  a3  +  h6 


(3-17) 


(3-18) 


where 


(1-s2) 

(1-t2) 


h$  and  h6  are  the  incompatible  Interpolation  functions. 
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ISOPARAMETRIC  ELEMENT 
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For  two-dimensional  analysis  the  strain-displacement 


equations  are 


EXX  "  3x  *  1  V k  uv<  +  h5.x  “l  +  h 


1,x  Ux1  "5.x  |  *  "6,x  “2 


du 

c  -  ~y . 


■  1  h1.y  Uy1  +  h5.y  “3  *  h6.y  “4 


Sux  .  3u, 


ay  *  lx  *  r  h(,y  “xl  +  E  h1.x  “y(  ♦  Hityai  +  h6>yo2  ♦  h5 >x„3  a 


\**k 


(3-19) 


Or  Equation  3-19  can  be  written  in  matrix  form  as 


c  -  §  u 


H,x  0  Ux 

°  l 

H,y  h,x  iy 


(3-20) 


3  x  12  matrix.  The  sobma, rices  In  equation  3-20  are  gi„„  „y 


S'x  '  thl.x  h2,x  »3,x  h4,x.h5.x-  hS,xJ 
aV  '  thl.y  h2.y  h3,y  h4>y'  h5,y.  h4.J 


(3-21) 


The  element  stiffness  matrix  is  given  by  the  following 


equation: 


-  "  «(  §T  5  B  dv 


(3-22) 


R-7215-2299 


A 

where  C  is  the  stress-strain  matrix.  The  above  equation  Is  Integrated 
numerical ly 

-  -  23  b|  C  B,  (3-23) 

This  stiffness  matrix  which  is  12  x  12  is  reduced  to  8  x  8  by  elimination  of 
the  four  incompatible  modes  before  assembling  in  the  global  stiffness  matrix. 

3.3.4  THrtEE-D I  MENS  I ONAL  ELEMENT 

For  an  arbitrary  eight-point  brick  element  shown  In  Figure  3*3,  the 
appropriate  displacement  approximations  are 

8 

“*  "  *  h9  °xl  * 

8 

uy  '  “yl  th9ayl  *  *10  V  "" 

8 

ui  *  *  u2l  +  h9  “zl  +  h10  az2  +  *11 

*  1«i 


°x3 

Oy3  (3-24) 

°z3 
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where 

h,  ■  1/8  (1  ♦  5)  (l  ♦  n)  (1  +  ;) 

li2  ■  1/8  (1  -  O  (1  ♦  n)  (1  +  c) 

1>3  -  1/8  (1  -  e)  <1  -  n)  (1  +  ?) 

h4  -  l/e  (1  ♦  t)  (1  -  n)  (1  ♦  0 

h5  •  1/8  (1  ♦  f)  (1  +  n)  (1  -  c)  (3-25) 

h6  ■  1/8  (1  -  {)  (1  ♦  n)  (1  -  0 

1*7  ■  1/8  (1  -  0  (1  -  n)  (1  -  0 

1*8  ■  1/8  (1  ♦  S)  (1  -  n)  (1  -  c) 

1*9  ■  (1  -  C2) 
h10  ■  (1  -  n2) 
h-,1  •  (1  -  c2) 

The  first  eight  are  the  standard  compatible  interpolation  functions.  The  last 
three  are  incompatible  and  are  associated  with  linear  shear  and  normal  strains. 
The  nine  incompatible  modes  are  eliminated  at  the  element  stiffness  level  by 
static  condensation. 


3.3.5  THICK  SHELL  ELEHENTS 

The  thick  shell  element  described  here  was  initially  developed  by 
Wilson,  et  al.,  Reference  3~3. 

This  shell  element  is  a  lb-node  curved  solid  element  shown  in 
Figure  3*4.  Each  node  has  three  unknown  displacements.  Therefore,  if  the 
shell  is  considered  as  a  two-dimensional  surface  there  are  six  unknowns  per 
point.  It  is  apparent  that  this  type  of  formulation  avoids  the  problems 
associated  with  the  sixth  degree  of  freedom--the  normal  rotation  is  set  to 
zero  when  certain  finite  elements  are  used  in  the  idealization  of  shells. 
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The  locations  of  the  nodes  are  defined  by  the  orthogonal,  right- 
handed  coordinate  system  (x,  y,  z)  which  is  referred  to  as  a  global  system. 

Within  the  element  a  local  coordinate  system  (£,  n»  O  has  been  chosen  such 
that  C.  n,  ;  vary  from  -1  to  +1;  (0,  0,  0)  Is  located  at  the  centroid  of  the 
element. 

The  local  and  global  coordinate  systems  are  related  through  a  set 
of  Interpolating  functions: 

16 

x  ■  E  h4  x4 
1-1  1  1 

16 

y  ■  E  h,  y,  (3-26) 

1-1  1  1 

16 

2  -  I  h4  zA 
1-1  1  1 

i 

/ 


Ik 


where 
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h1  ■  V8  (1  ♦«)  <1  +  n)  (1  ♦  c)  <{  ♦  n  -  1) 

f2  •  1/8  (l  -  £)  (1  +  n)  (1  +  ;)(-{  ♦  i,  -  1) 

h3  *  V8  ('  -  0  (1  -  n)  (1  ♦  ;)(-£  -  n  -  l) 

h4  ■  VS  (1  +  «)  (1  -  n)  (1  ♦  c)  (£  -  „  .  i) 

h5  ■  V8  (1  ♦  ?)  0  ♦  n)  (1  -  t)  (£  ♦  „  .  ,) 

h6  •  1/8  (i  -  £)  (1  ♦  n)  (1  -  ;)(-£  +  n  -  l) 
h7  *  '/8  (1  -  «:  (1  -  n)  (1  -  ;)(.£  .  n  .  i)  (3'2; 

ha  ■  !/s  (1  +  0  (1  -  n)  <1  -  c)  (e  -  n  -  i) 

"s  ■  1/4  (1  -  £2)  0  +  n  )  (l  +  c) 

hio  "  '/4  (i  -  t  )  (1  -  n2)  (l  +  c) 

h„  ■  1/4  (i  -  £2)  (1  .  n  )  (1  +  c) 

h12  *  1/4  (1  ♦  £  )  (1  -  n2)  (1  +  £) 

h,3  •  1/4  (1  -  £2)  (1  ♦  n  )  (1  -  £) 

8.4  ■  1/4  (1  -  £  )  (1  -  n2)  (,  .  C) 

8.5  *  1/4  (1  -  £2)  (1  -  n  )  (1  -  0 

8.6  -  1/4  (1.4  £  )  {,  -  r,2)  (1  -  C) 

Th.  displacements  within  the  element  are  assumed  to  be  of  the  follow 

j  form: 

16 

'  iJ1  h1  Ux1  +  h17  “xl  +  h18  ax2  +  h19  ax3  +  h20  ax4  +  h21  °x5 
1C 

"i  '  h1  Uyi  +  h17  ayl  +  h18  °y2  +  h19  ay3  +  h20  ay4  +  h21  ay5  (3-28) 
16 

”,  "  iJ1  h1  uzi  +  h17  az1  +  h18  az2  +  h19  az3  +  h20  az4  +  h21  rz5 
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whe  re 

h17  ■  E  0  -  E2) 
^18  *  11  ^  '  n2) 


^20  *  En  (1  -  £2) 
h2l  *  nc  (1  -  n2) 

The  motivation  for  addition  of  the  Interpolation  functions  h  to 
h2J  Is  to  Increase  the  capability  of  the  element  In  producing  closer  approx! 
(nations  to  the  exact  displacements  under  simple  loadings,  thereby  increasing 
the  convergence  to  exact  solution.  The  Incompatible  Interpolation  functions 

h| 7  t0  h21  have  zero  va,ues  at  the  n°des  and  produce  Incompatibilities  in 
displacement  field  along  the  intereiement  boundaries. 

3. A  JOINT  FINITE  ELEMENT 


The  joint  element  is  intended  to  represent  the  rock  joints,  faults, 
interfaces  and  similar  discontinuities  in  continuum  systems.  The  joint 
element  has  the  capability  of  representing  the  main  characteristics  of  the 
deformation  behavior  of  the  rock  joints  such  as  debonding  and  slip.  The 
term  debonding  means  the  ability  of  separation  of  the  two  blocks  of  continuum 
adjacent  to  the  joint  surface  which  were  initially  in, contact.  Subsequent 
contact  can  also  develop  by  the  movement  of  the  two  blocks  towards  each  other. 
The  term  slip  means  the. relative  motion  along  the  joint  surface  or  fault  when 
the  shearing  force  exceeds  the  shear  strength  of  the  joint. 

Previous  attempts  have  been  made  to  develop  discrete  elements  to 
represent  the  joint  behavior.  Goodman,  Taylor  and  Brekke  (Reference  3-2) 
developed  a  simple  rectangular,  two-dimensional  element  with  «*ight  degrees 
of  freedom.  This  element  has  no  thickness,  and  therefore  the  adjacent  blocks 
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of  continuum  elements  can  penetrate  into  each  other.  Zienkiewicz,  et  a)., 
(Reference  3”^)  advocate  the  use  of  cont inuum  isoparametric  elements  with  a 
simple  nonlinear  material  property  for  shear  and  normal  stresses,  assuming 
uniform  strain  in  the  thickness  direction.  Numerical  difficulties  may  arise 
from  ill  conditioning  of  the  stiffness  matrix  due  to  very  large  off-diagonal 
terms  or  very  small  diagonal  terms  which  are  generated  by  these  elements  in 
certain  cases. 

To  avoid  such  numerical  problems  a  new  joint  element  is  developed 
below,  which  uses  relative  displacements  as  the  independent  degrees  of  freedom. 
The  displacement  degrees  of  freedom  of  one  side  of  the  slip  surface  are 
transformed  into  the  relative  displacements  between  the  two  sides  of  the  slip 
surface.  The  transformation  relations  are  as  follows: 


T 

u  . 


u 


u 


u 


x  i 
T 

yi 

T 

xj 

T 

yj 


uxl  +  4uxl 


u  ,  +  Au  . 
yi  yi 


u  ,  +  Au  . 
xk  xj 


u  ,  +  Au  . 

yk  yj 


The  superscripts  T  and  B  refer  to  the  top  and  bottom  elements 
with  respect  to  the  slip  surface  respectively.  As  shown  in  Figure  3-5,  those 
degrees  of  freedom  of  the  upper  element  which  are  on  the  slip  surface  are 
transformed  but  the  degrees  of  freedom  of  the  lower  element  are  the  original 
displacement  quantities. 
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UPPER 
CONTI 
ELEMENT 


* 

JOINT  ELEMENT 

\ 


LOWER  CONTINUUM 
ELEMENT 


■ - -«►  x 

(GLOBAL  COORDINATES) 


xi 

uT. 

y 

T 

j  . 

*j 

T 

Jyj 


u  .  +  Au  . 
xl  xl 


uyl  +  4uyi 


u  ,  +  Au  . 
xk  xj 


u  .  +  Au  . 
yk  yj 


FIGURE  3-5.  GEOMETRY  OF  JOINT  ELEMENT 
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The  joint  element  is  assumed  to  have  the  relative  displacements  as 
the  degrees  of  freedom.  For  example.  In  a  two-dimensional  problem  the  joint 
element  will  have  four  degrees  of  freedom  (Figure  3-6).  The  relative  normal 

and  tangential  displacements,  4up  and  4us,  are  assumed  to  vary  linearly  along 
the  element  as  follows 


Au  =  h  Au  +  h.Au  . 
n  f  ni  J  nj 


Au  =  h.Au  .  +  h.Au  . 
s  i  si  j  sj 


(3-30) 


where  the  h.  and  h.  are  the  I 


inear  interpolation  functions 


h,  -  i(i  -  e) 

hj  ■  l('  +  «) 


(3-31) 


and  iun.,  iu^j,  4us|  and  iu^  are  the  nodal  point  values  of  the  relati\ 
sp lacements .  The  joint  element  is  assumed  to  have  only  two  strain  com¬ 
ponents;  en  .  normal  strain,  and  es  =  shear  strain.  These  two  strain 

exponents  are  related  to  the  relative  displacements  through  the  following 
relations. 


e  =  -  Au 

n  t  n 


es  "  TAus 


(3-32) 
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The  substitution  of  Equations  3“30  and  3-31  into  Equation  3-32  results 
in  the  strain-displacement  relation  for  the  element 


(3-33) 


The  stresses  and  the  strains  are  related  through  the  following 
material  property  matrix  C. 


g  ■  C  e 


(3-34) 


In  general  the  above  suess/stra i n  relationship  for  rock  joints  is 
nonlinear,  the  details  of  which  are  given  in  Section  4. 

The  stiffness  matrix  for  the  joint  element  is  formed  in  n-s 
coordinate  system; 


•ns 


f  BTC  B 
-  -  - 


dv 


(3-35) 
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and  transformed  to  the  x-y  coordinate  system  as  follows 

k  -  TTk  T 
-  -ns  - 

where  J  is  the  transformation  matrix  containing  the  direction  cosines. 


(3-36) 


2 (A i  "  2V  2(A3  +  B2)  (A,  -  2B , )  (A3  +  B2) 

2(A2  +  2B j )  (A3  +  B2)  (A2  +  2Bj) 

c  ♦  •  2(A|  "  2V  2 (Ao  +  BJ 

Symmetric  11  32 


2(A2  +  2B2) 


(3-37) 


where 


0  9 

Css=-  * 
c  b2+C  a2 

ss  nn 


(C  -  C  )  ab 
nn  ss' 


f.  ab 
ns 


c„5U2  -  b2> 


a  '  r(xj  *  xi> 
"  '  f  (*j  ‘  »t> 
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SECTION 

REPRESENTATION  OF  PROPERTIES  OF  ROCK,  INCLUDING  ANISOTROPY, 
INELASTICITY,  RATE  EFFECTS  AND 
PROPERTIES  OF  FAULTS  OR  JOINTS 

The  first  part  of  this  section  describes  homogeneous  properties  of 
rock  which  are  available  in  the  AA  computer  program.  As  used  here,  "homogeneous1' 
refers  to  properties  which  can  reasonably  be  averaged  over  several  feet,  which 
is  a  typical  dimension  of  a  finite  element  in  applications  to  mining  engineer¬ 
ing.  Inhomogeneous  properties  of  rock  masses,  such  as  those  caused  by  faulting, 
arc  treated  by  a  separate  procedure  which  is  described  in  the  second  part. 

The  topics  which  are  covered  beiow  include  the  following: 

a.  Inelasticity 

1 .  Variable  modulus 

2.  Variable  modulus  with  perfect  plasticity 

3.  Variable  modulus  with  perfectly  plastic  fracture 
criterion  and  strain  hardening  cap 

b.  Anisotropy 

1.  Variable  modulus  with  anisotropic  fracture  criterion 
based  on  the  hypothesis  of  Jaeger  (plane  geometry  only) 

2.  Variable  modulus  wi th  anisotropic  yield  criterion  based 
on  the  hypothesis  of  Niii 

c.  Rate  Effects  (Isotropic  Oniy) 

1.  Creep  (series  of  Kelvin  elements) 

2.  Viscoplasticity  (based  on  work  of  Perzyna) 


3* 


Joint  Properties 
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d. 


'■  Dilatant 
2.  Nondi i I tant 


*»•  I  homogeneous  properties 

T'lese  models,  which  are  also  summarized  in  Table  <♦-!,  afford  some 
material  description  which  probably  cannot  be  affectively  used  at  present 
owlnp  to  the  leek  of  experimental  data  on  rock  from  particular  mines.  As 
analysis  Increases  In  effectiveness.  It  Is  expected  that  more  data  will  be 
routinely  obtained,  thus  enabling  more  sophisticated  material  models  to  be 
used.  Advantages  and  disadvantages  of  these  models  are  summarized  In 
Tab.e  4-2.  Typical  laboratory  properties  for  rocks  are  shown  In  Figure  4-1. 


AuLE  4-1.  SUMMARY  OF  AVAILABLE  MATERIAL  PROPERTIES 


Materia  1 

Characteristics 

Type  of  Model 

/ 

/• 

/ 

S* 

V* 

/ 

/ 

gN 

/ 

Constant  elastic  moduli 

X 

0 

X 

0 

Variable  modul 1 

X 

0 

X,0 

Elast lc/plast Ic  (fixed 
static  yield  surface) 

Constant  moduli 

X 

0 

X 

0 

Variable  moduli 

X 

0 

X 

0 

Cap  model 

X 

X 

X  options  or  0  options  may  be  used  simultaneously,  but 
x  and  0  options  may  not  be  mixed,  except  where  noted. 
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TABLE  4-2.  ADVANTAGES  AND  DISADVANTAGES  OF  EACH  MODEL 


Advantages  Disadvantages 

A.  Elastic-Ideal ly  Plastic 


Simple  to  f I t  to  data 
Approximates  most  features  of  data 
G  -  Const.  G  ■  G(P  )  and  associated 
flow  rule  theoretically  correct 


May  not  fit  all  available  data 
Cannot  match  triaxlai  test 
Other  treatments  of  G  can  lead 
to  possible  paths  of  energy 
generation 

For  nonassoclated  flow  rule  no 
general  uniqueness  theorem 


B.  Variable  Moduli 


Best  fit  of  data  Restricted  to  near-proportional 

Only  model  with  repeated  hysteresis  loading  (in  shear) 


within  failure  envelope 

Ideal  for  finite  element 
Computationally  simple 
Relatively  easy  to  fit 

For  nonproport  Ion  a  1  loading  paths 
no  uniqueness  theorem 

Additional  quantity  must  be  stored 
at  each  grid  point 

C. 

Cap  Model 

Satisfies  ail  rigorous  theoretical 
requi remen ts 

Reasonably  good  fit  of  data 
Effective  control  of  dliatancy 

Indirect  approach  needed  to  fit 
data 

Relatively  complicated 

Additional  quantity,  the  strain 
hardening  parameter,  must  be 
stored  at  each  grid  point 

D. 

VI scoelastic 

Simple  to  fit  to  data 
Approximates  features 
of  data  for  some  rocks 

Requires  sophisticated  testing 
to  define  viscous  coefficients 
for  mut ti-axiai  loading. 

Does  not  account  for  deterioration  in 
strength  with  time 

E. 

Viscoplastic 

Approximates  some  features  of  Requires  sophisticated  testing 

data  Including  shear  strength, 
stick-slip  phenomenon 
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TRI AXIAL  COMPRESSION 


Ml 

TRI AXIAL  COMPRESSION 


AV/V 

AJAJ7S1 


FIGURE  4-1.  TYPICAL  LABORATORY  DATA  ON  ROCK  FROM  WHICH  CONSTITUTIVE 
EQUATIONS  ARE  OERIVEO 
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The  programming  of  the  material  property  subroutines  Is  arranged  to 
provide  maximum  flexibility  and  ease  In  changing  properties.  This  Is  done  by 
performing  each  separate  function  In  the  material  properties  description  by  a 
separate  subroutine.  Thus  e»arate  subroutines  are  provided  for  the  following 
purposes : 


Computing  variable  moduli 

Computing  derivatives  of  the  yield  functions  with  respect 
to  its  arguments. 

Testing  for  yielding  or  fracture 

Adjusting  stresses  for  viscoelastic  relaxation. 

Adjusting  stresses  for  viscoplastic  relaxation. 

Transforming  strain  Increments  to  principal  axes  of  orthotrophy 
and  transforming  the  matrix  of  stress/straln  coefficients  to 
the  global  axes. 

There  is  considerable  interdependence  between  the  inelastic  and  anisotropic 
capabilities.  V/herever  possible  this  Interdependence  is  used  to  economize  on 
programming.  Logic  diagrams  for  each  subroutine  In  the  homogeneous  material 
property  package,  are  given  In  Appendix  A. 


*•!•!  IHELASTI CITY  FOR  ISOTROPIC  MATERIALS— VARIABLE  MODULUS 

Inelasticity  In  Isotropic  materials  Is  represented  through  variable 
bulk  and  shear  moduli  and  throi-jh  plasticity  theory.  The  bulk  modulus  B 
Is  assumed  to  depend  on  the  current  value  of  elastic  volumetric  strain  p  and 
its  previous  maximum  value  u 

max 


FOR  LOADING  (0  >  u  £  u^) 


B 


(U  -  B  )  exp 
mo 
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FOR  UNLOADING  OR  RELOADING  (0  <  u  < 

B  '  Bu  *  (Bm  ‘  Bu)(5j)  (*-*) 

where 


(Bo  +  (Bm  •  Bo> 

Bu  ■  the  lesser  of  v 

IB 

V  m 

FOR  LOADING  OR  UHLOADING/RELOADING  IN  TENSION  (u  <  0) 

B  -  Bt  (1,-3) 

Application  of  this  model  to  a  granitic  rock  is  Illustrated  in 
Figure  k-2.  Specific  parameters  for  this  rock  are 

Bm  -  7-6  x  I06  psi 

Bq  -  1.205xl06psl 
u,  -  0.0275 
u2  -  0.05 

The  shear  modulus  G  is  also  assumed  to  depend  on  u  and  u 

max 

FOR  LOADING  (0  <  u  <  u) 

*  WX 

G  ■  Gm  •  (Gm  •  t0)  “O  (j-j)  <»•« 
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(a)  MODEL  BULK  MODULUS  FOR  LOADING 
AND  UNLOADING 


A/ 

\A//  .  . 

T  O.Ol  0,02  0.03  0.04  0.05  0.06  0.07  0.08 
P 

AJAI051 

(b)  MODEL  HYDROSTAT  FOR  LOADING 
AND  UNLOADING 


I  -  - 


30.00 


25.00  — f 


••  / 


J 

/-/ 


*  M 

■  *J  L***  ml 

07  ,  **  L*  *  .  -  r  — 


*  .  PRESENT  £ 
£_J  MODEL 


/ 

15.00  I  *  / 

h 

10.00  SJ / 


-  UNIAXIAL 

STRAIN 

-  HYDROSTAT 


•  //  :ss 

\£a  I  1  1  1  ]n 

0  0.02  0.04  0.06  0.08  0.10  0.12  0 


0  0.2  0.4  0.6  0,8  1.0  0  0.02  0.04  0.06  0.08  0.10  0.12 

PRESSURE,  Kb  V 

(c)  LOW  PRESSURE  BULK  MODULUS  COMPARED  (d)  MODEL  HUGONIOT  AND  HYDROSTAT 
WITH  DATA  (REFERENCE  4-1)  COMPARED  WITH  DATA 


FIGURE  4-2.  MODEL  BULK  MODULUS,  HYDROSTAT  AND  YIELD  POINT  IN  UNIAXIAL 
STRAIN  COMPARED  WITH  DATA 


FOR  UNLOADING/RELOADING  (0  <  \i  <  u  ) 

"  'fnax' 
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G  =  Maximum  previous  value  of  G 


(4-5) 


FOR  LOADING  OR  UNLOADING/RELOADING  IN  TENSION  (u  <  0) 


G  =  G. 


(4-6) 


Application  of  this  model  to  cracked  and  uncracked  granitic  rocks 
is  illustrated  in  Figure  4-3.  Specific  parameters  for  these  rocks  are: 

Gm  =  4.35  x  I06  psi 

G  =  0 

o 

u3  =  0.005 


n  The  incremental  stress  da.j  is  related  to  the  increments  com¬ 

ponent  of  elastic  strain  de?^  by  the  following  expression: 


where 


.  dcu  r  (B’fG)  dekk  6ij +  2G(deTj) 


6ij  =  0  if  1  +  Ji  "  1  'f  '■  -  J' 


(4  /) 


Theoretical  guidance  on  the  appropriate  functions  for  B  and  G 
is  provided  by  Walsh  (References  4-7.  4-8),  who  postulates  that  the  effective 
modulus  differs  from  the  intrinsic  modulus  due  to  cracks  and  pores.  As  these 
are  closed  by  increasing  pressure,  the  effective  modulus  tends  toward  the 
consolidated  value.  Walsh's  work  contains  parameters  which  are  not  retained 


r-7215-2299 


PRESSURE,  Kb 


FIGURE  4-3.  SHEAR  HODULUS  VERSUS  PRESSURE  FOR  NTS  GRANITE 
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in  the  following  empirical  expressions  for  the  effective  bulk  modulus.  How¬ 
ever,  the  basic  concept  is  retained.  Also,  the  present  model  for  the  effective 
shear  modulus  merely  follows  Walsh's  concept.  The  idea  of  coupling  the  shear 
stiffness  to  volumetric  strain  is  proposed  in  Reference  4-9  and  carries  the 
danger  that  energy  might  be  extracted  from  the  model  by  hydrostatic  com¬ 
pression,  followed  by  shearing,  followed  by  releasing  the  pressure  and  finally 
by  releasing  the  shear.  This  danger  is  avoided  by  assuming  that  friction 
prevents  cracks  from  reopening  during  unloading  so  that  the  largest  value  of 
G  reached  on  loading  is  retained  during  subsequent  unloading/reloading. 

Under  these  conditions  a  material  may  dissipate  energy  in  shear  during 
loading  and  unloading  cycles  but  can  never  produce  additional  energy. 

4.1.2  INELASTICITY  FOR  ISOTROPIC  MATERIALS— VARIABLE  MODULI  WITH  PLASTICITY 

The  present  adaptation  of  plasticity  theory  is  based  on  work  of 
References  4-9  through  4-13.  The  model  consists  of  a  yield  criterion 

f(a.j,  L)  =  0  (4-8) 

where  L  is  a  function  of  plastic  strain, 

and  a  plastic  flow  rule  in  which  f  is  regarded  as  a  potential  function 


d 


3f 


3a. 


U 


(4-9) 


The  incremental  stress  is  related  to  the  elastic  component  of  incremental 
strain  by  Equation  4-7.  Defining 


de?.  =  dc. .  -  dc? . 
•J  'J  ij 


(4-10) 
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Substituting  Equations  4-9  and  4-10  into  Equation  4-7  leads  to 

doij  ■  AKk-Ai^)sij +  4cu  -Ai%) 

where 

A  =  li  -  jG 

If  the  yield  criterion  is  satisfied,  the  stress  state  must  lie  on  the  surface 
defined  by  f  in  Equation  4-8.  The  mathematical  statement  of  this  constraint 
i  s 


af  ,  af  .. 

df  «  - -  do.  .  +  —  dL 

3o. .  ij  aL 
■J 


Substituting  Equation  4-11  into  4-12  permits  solutions  for  A 


X^kk1  f«  *  2G  dei  if  ii 
>fkkfM  *  2G  fijfij  +  » 


(4-12) 


(4-13) 


where  R  is  a  strain-hardening  function  to  be  defined  below. 

Substitution  of  Equation  4-13  into  4-11  expresses  the  stress  increment  in 
terms  of  the  strain  increment. 


A 
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Specific  functional  forms  have  been  assumed  for  f.  These  contain 
empirical  constants  whose  values  can  be  selected  to  match  data  for  a  specific 
material.  The  forms  are 


Polynomial  1  In  o, 


V7!  ‘  £,  v 


V7!- 


0  J,  >  b 


0  J,  <  b 


(A-lA) 


Polynomial  2  in  o. . 


fi  <0i  j5 


V3!  -  a.  '  -  ('  •  4) 

V7!  ■  (a!  +  a2} 


0  J 1  >  b 


0  J,  1  b 


(*-15) 


Cap  (to  be  used  with  Polynomial  2) 


f2  -  (J,  -  V)2  *  P2(j-  -  a)  -  0 


(*-16) 


in  which 


V  -  L  +  P  X(L)  X  * (L) 


(X(L)]2{i  +  P2[X'(L)]2| 


(*-17) 


(*-18) 


-('-lr)2l 


A  +  C 


L  S  B 


L  >  13 


(*-19a) 


*5 
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A 

ir) 

L  >  B 

The  hardening  parameter,  L,  is 


where 


L  - 


iP2dt 


g 


w  1 


2 


(4-19b) 


(4-20) 


(4-21) 


(4-22) 


The  hardening  parameter  R  is 


(4-23) 


where 

°1’  °2*  °3  are  PriPc{PaI  stresses 

and 

R  -  0  if  f  -  f,. 

The  cap  parameters  and  stress/strain  relations  produced  by  the  cap 
model  are  shown  in  Figures  4-4  through  4-6.  Data  on  strength  for  granite 
containing  various  degrees  of  cracking  are  shown  in  Figure  4-7,  which  illus¬ 
trate  the  adequacy  of  the  assumed  fracture  criteria. 
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RADIAL  STRESS,  KSI 


FIGURE  4-5. 
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CAP  MODEL  FIT  FOR 
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-  VOLUMETRIC  STRAIN 

-  SHEAR  STRAIN 


o3  -  250  KSI 


(SEE  NOTE  ON 
ATTACHED  SHEET) 


STRAIN,  PERCENT 


FIGURE  4-6.  TRIAXIAL  STRESS  BEHAVIOR,  CHAMBER  PRESSURE  -  0,  50,  100  200  KSI 
CAP  MODEL  FIT  FOR  GRANITIC  MATERIAL 


FIGURE  J»-7.  YIELD  STRENGTH  VERSUS  MEAN  STRESS  FOR  NTS  GRANITE 
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The  Incremental  stress/strain  equations  are  expressed  In  matrix  form 
as  follows: 

{do}  ■  (C]{dc}  (4-24) 

wnere  dc  Is  the  total  increment  of  strain.  The  C  matrix  thus  contains 
generalized  tangent  moduli  and  can  be  used  In  forming  the  element  stiffness 
matrices.  For  the  models  descr'bed  above,  the  C  matrix  is  as  follows: 
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The  C 


matrix  Is  clearly  composed  of  elastic  and  Inelastic 


parts 


■e  -  r  P 


(4-26) 


The  C*  and  Cp  matrices  are  competed  separately.  The  reason  for  this  Is 

efficiency  In  tr.a,lng  Isotropic  .„d  anisotropic  materials  with  ,h.  same 
Fortran  statements. 

r.fo  ,  !"  “T'V  n°"lin'i,r  pr°bl""5  11  15  *>»“»«  to  avoid  time-consuming 

ormulatlon  of  the  stiffness  matrix  by  Introducing  nonlinearity  through 

t  e  oad  vector.  The  method  is  an  extension  of  the  egul  librium  equations 
given  In  R.f.r,„«  1,-29.  In  the  following  equation,  time  Is  used  as  a 
parameter.  Number  of  load  step  couid  be  used  Instead. 

At  time  t,  -he  total  change  in  complimentary  strain  energy  is 
equal  to  the  change  in  complimentary  work  done  by  nodal  point  forces. 


T-t  f 

£  / 

T"0  JMO I 


<C>T  {do)T  dV 


£  <u>  (dP) 

T-0  T  T 


where 


<=>,  (do) 
<u>,  (dP) 
V 


Element  strain  and  stress  Increment 
Nodal  displacement  and  force  vector 
Volume  of  finite  elvnent 
Arbitrary  Instant  of  time 


The  strain/dl spiacement  relation  Is 


(e)  -  (Bj{u) 

-  <u>[B]T 


(4-27) 


(4-28) 
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The  stress  increment  in  an  elastic/plastic  material  may  be  expressed  by 
rewriting  Equation  4-24  as  follows. 


{da}  ■  [C] { de>  -  (dap} 


(4-29) 


elastic  correction 


where  now  C  •  C 


Substituting  Equations  2-28  and  2-29  into  Equation  2-27, 


£  f  <u>T  [B]T{[C] [B]{du>x  -  (dop)T|dV  -  S  <u>T  (dP} 
T-0  •'vol  T  T~° 

Noting  that  <U>T  maY  be  e 1 i m i na ted  from  both  sides  and  that 

f  [B]T[C] [B]dV  -  [K] 

•'vol 


(4-30) 


(4-31) 


where  [K]  is  the  elastic  stiffness  matrix,  Equation  2-30  may  be  rewritten  as 


T(ik,{^  ■  L IB,T,d^ dV) 


(4-32) 


+  [K](du}  -  f  [B]T{da  }  dV  -  £  {P)T 

Jvo i  V  At  t=0 

where  {duk  ,  jda  |  equal  change  in  u  and  o  during  the  interval 
At  «  p'At  K 

t  -  At  to  t. 

Performing  the  indicated  summation,  assuming  stress  to  be  constant 
throughout  the  element  and  defining 


E  -  £  [KHdu) 

T  T-0  T 


(4-33) 
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[K] {du} 

At 


(P,t  -  +  t“lTK}  V 

'  H,t-At 


(4-34) 


For  small  increments  of  stress  and  time  in  a  step-by-step  integration,  the 
second  term  on  the  left  hand  side  is  neglected. 


The  expression  which  is  used  for  {E}  in  the  computer  program  is 
derived  as  follows.  The  recoverable  work  done  on  an  element  is  equated  to 
the  elastic  strain  energy  stored  in  the  element  by  the  following  equation. 


l<u>{E)  =  -<£>[C]{e}V 

or 

(E)  =  [BT][C]{E}V 

Thus  Equation  4-34  may  be  rewritten 

[K](du)it  .  ff)t  -  [BTj |[C]{eJt.it  -{cp}  lV 
Following  Equation  4-29 

d0ij  ■  (XdEkkSij  ♦  2G  dEij)  ‘  (x  5^  S!j  +  2G  3?fj)  A 


doij  ■  xdEkksij  +  2G  dEij  -  dEij 


(4-35) 


(4-36) 


(4-37) 


(4-38) 


(4-39) 
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where 


Xde 


9f 


kk  9 o 


+  2G  d~. 


9f 


U 


i  j  9  o .  . 

'JL 


x  -3JL  _if_  +  ?G  JlL 

Sokk  9ou  90  -  90 - 


•  j 


•  j 


o 


if  f  -  o 


(4-4o) 


if  f  <  0 


After  each  integration  step, 


Au. 

i 


is  given; 


de. .  is  found . 
•J 


Based  on  de. ^  and  stress  at  previous  time,  f  is  checked 
and  da.j  is  calculated  according  to  whether  material  is 
elastic  or  plastic  by  adding  contribution  for  elastic  part 
to  contribution  from  plastic  part. 


Stiffness  is  always  based  on  elastic,  parameters,  i.e.,  X,  G. 
Plasticity  is  introduced  through  updating  of  stress  increment. 
Hence,  there  is  need  to  update  stiffness  matrix. 


4.1.3  INELASTICITY  FOR  ANISOTROPIC  MATERIALS— THEORY  OF  JAEGER 

The  fracture  of  anisotropic  rocks  is  the  subject  of  several  failure 
theories.  The  Walsh-Brace  theory  (Reference  4-14)  assumes  that  failure  is 
tensile  in  nature  and  that  it  is  influenced  by  the  presence  of  preexisting 
cracks.  Some  of  the  cracks  are  assumed  to  be  small  and  randomly  oriented 
while  others  are  long  and  have  preferred  directions.  Extension  of  the  cracks 
is  postulated  to  occur  when  the  Griffith  criterion  (Reference  4-15),  as 
nodified  by  McCiintcck  and  Walsh  (Reference  4-16)  to  account  for  friction  on 
the  crack  faces, is  satisfied. 


I  n 

shear  either 
according  to 


contrast,  Jaeger  (Reference  4-17)  assumes  the  material  to  fail  in 
along  a  single  plane  of  weakness  or  within  the  matrix  material 
a  Mohr-Coul omb  type  of  failure  criterion  of  the  form 


t  =  a  oa, 

o  1 


(4-38) 
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Normal  stress  on  plane  of  fracture 
Cohesion,  angle  of  friction 
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where 


a 

3  >  3 


The  theory  is  expressed  in  terms  of  principal  stresses  by  an  axis  rotation 
as  follows 


2ao  ~  2al°3 

ai  •  Vaf  - 1 


(**-39) 


Improved  agreement  with  experiment  is  obtained  If  aQ  is  assumed  to  vary 

ao  =  a2  "  a3  cos  (2(*  '  *))  (**-**0) 

where 

6  =  Counterclockwise  angle  from  the  direction  of  the  major  prin¬ 

cipal  stress  (cfj)  to  the  direction  of  the  bedding  planes 

5  -  The  orientation  of  3  at  which  a  is  minimum.  Usually 

assumed  to  be  equal  to  30  deg  ° 

The  angle  3  is  shown  in  Figure  **-8  along  with  angles  relating  it  to  global 
directions  in  the  finite  element  formulation.  McLamore  and  Grey  (Refer¬ 
ence  **-18)  have  obtained  satisfactory  agreement  between  experimental 
data  on  strengths  of  slates  and  shales  using  a  modification  of  Equation  **-25 
as  follows: 


aQ  ■  a2  -  (cos  2(5 

ao  =  a4  "  a5(cos  2(5 

Some  of  their  results  and  those  of 
Figure  **-9. 


-  e)]n 

for 

0  <  5  <  3 

-  e)]n 

for 

3  <  5  <  90' 

Brace  and  Walsh  are  illustrated  in 


(**-**1) 
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X  -  Y  -  Z  =  GLOBAL  AXES 
1-2-3  -  PRINCIPAL  STRESS  AXES 


FIGURE  i*-8.  ORIENTATION  OF  PLANES  OF  WEAKNESS  DEFINING  ANISOTROPIC  BEHAVIOR 
ACCOROING  TO  JAEGER 
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FIGURE  4-9.  COMPARISON  BETWEEN  FAILURE  THEORIES  AND  EXPERIMENTAL  DATA 
FOR  ANISOTROPIC  ROCK  (REFERENCE  2-18) 
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The  present  application  of  the  hypotheses  in  Reference  4-18  is 
limited  to  plane  geometry.  It  does  not  take  into  account  the  effect  of  the 
intermediate  principal  stress  which  is  sho-'n  in  Reference  4-19  to  play  an 
important  role  in  fracture  of  some  types  of  rocks.  More  work  is  needed  to 
remove  these  restrictions. 

The  first  step  is  to  determine  the  magnitude  of  the  principal 
stresses  o^,  and  their  direction  as  specified  by  0,  the  counter¬ 
clockwise  angle  from  the  +X  global  axis. 


(4-42) 


(4-431 


Bilinear  stress/strain  relations  are  specified  by  the  user  in  terms 
of  Young's  moduli  and  Poisson's  ratio  in  directions  parallel  and  perpendicular 
to  the  bedding  planes  (oriented  at  a  relative  to  global  axes).  Thus, 
experimental  data  is  required  from  specimens  cut  orthogonal  to  bedding  plane 
and  at  angles  other  than  90°.  The  computer  program  transforms  the  various 
E  and  v  to  the  principal  directions  of  stress  and  modifies  them  to  account 
for  fracture.  These  parameters,  Ep  v 12>  v^,  Ep  Vp,  ,  etc.,  are 
assembled  into  a  matrix  relating  incremental  stress  and  strain  in  principal 
stress  axes.  The  relationship  between  incremental  stress  and  incremental 
strain  expressed  in  the  principal  axes  of  anisotropy  (principal  stress  axes) 
is  shown  in  Equation  2-24  where  C  is  given  by  Equation  4-44.  The  matrix  is 
then  transformed  through  the  angle  0  into  global  coordinates  for  inclusion 
in  the  element  stiffness  matrix.  An  illustration  of  the  bilinear  Young's 
modulus  approach  is  superposed  on  data  in  Figure  4-10. 
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FIGURE  4-10. 


STRESS-STRAIN  CURVES  FOR  GREEN  RIVER  SHALE-2  FOR 
CONFINING  PRESSURES,  6-10  DEG  (REFERENCE  2-18) 


VARIOUS 


A 


£l(l  •  *23 V 


£3(VI3  +  *32*21  ) 
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E1(*I3  +  *,2*23) 
£2(*23  +  *21*13) 
£3°  '  *12*2)) 


(4-M) 


where 


0  “  1  -  v 


23  32  •  12*21  -  *13*3,  -  *1 2*23*3 ,  '  *13” 


H  31  13  21v32 

IN£USTICITV  AIIISOTROPIC  rock-theory  of  hill 

tne  stress  components  are  expressed  in  rh«  •  , 

axes  of  anisotropy  (no,  necessarily  principal  stress  axes)  ,H  ,7 
can  be  expressed  by  ’  the  yie  d  criterion 


where 


f(V3F-  Jf)  ■  0 


H-45) 


J2‘"’  '  ClKt  •  %„)2  *  C2(om  -  ou)2  ♦  C3(o,c  -  „  )2 


C*°«  +  Vn?  *  c6°« 


Ji'  =  c7a«  =  c8ann  +  c9a 
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In  which  it  is  assumed  that  the  C,  n,  E  axis  system  coincides  with  the  principal 
direction  of  anisotropy,  f  may  be  used  as  the  potential  function  described 
above.  Adaptation  of  this  theory  to  rock  is  described  in  Reference  4-20. 

The  elastic  behavior  of  the  material  may  be  prescribed  to  be  either 
isotropic  or  anisotropi  .  If  it  is  isotropic,  the  quantities  B  and  G  may 
be  used.  If  i  is  anisotropic.  Young's  moduli  and  Poisson  ratios  Ej,  v12, 
v^,  etc.  are  specified  in  the  principal  directions  of  anisotropy.  The  C 
matrix  (Equation  4-24)  which  relates  incremental  stress  to  incremental  strain 
is  thus  initially  expressed  in  the  principal  axes  of  anisotropy  and  is  sub¬ 
sequently  transformed  to  global  directions  of  the  finite  element  mesh. 

4.1.5  RATE  EFFECTS— VISCOPLASTICITY 

This  method  of  incorporating  rate  sensitivity  equations  is  based  on 
Perzyna's  elastic-viscoplastic  model  (Reference  4-21)  which  is  a  generaliza¬ 
tion  of  an  earlier  model  proposed  by  Hohenemser  and  Prager  (Reference  4-22). 

An  adaptation  of  the  cap  model  described  above  for  viscoplasticity  is  described 
in  Reference  4-23.  The  present  model  is  taken  from  Reference  4-24. 

A  linear  elastic,  rate  independent  region  is  bounded  by  a  static 
yield  cri ter  ion 

f(J,,  Jp  <  0  (4-46) 

within  which  Hooke's  Law  applies.  If  the  static  yield  criterion  is  satisfied 
or  exceeded 

f(J,,Jp  >  0  (4-47) 

A  viscoplastic  strain  rate  is  assumed  to  develop  according  to  the  following 
flow  rule. 


Y*(f> 

•J 


(4-48) 
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where 


4>(f)  -  A  function  of  the  static  yield  criterion  f 

F  -  Assumed  potential  function.  Presently,  a  nonassociated 

flow  rule  is  used  in  which  F  ■  and  differentiation 

is  performed  with  respect  to  the  stress  deviators  o' 

ij 

Y  “  Empirical  viscoplastic  parameter 


It  is  assumed  in  the  present  work  that 


♦  (f) 


FT 


\J2 


a  j;-1 

n  1 


n«l 


VT 


-  1 


J,  >  b 


-  1 


J,  s  b 


for 


V'T"  £  a  j"  1  >  0  J  >  b 

n- 1  1 


‘yp J  '  a5  >  0 


J 1  <  b 


(4-49) 


(4-50) 


If  f  <  0,  elastic  inviscid  stress/strain  relations  apply. 


Making  use  of  Houke's  law 


=  *=kkSij  +  2G  *Tj 


(4-51) 


and  expressing  the  elastic  deviatoric  strain  rate  by 


•  e ' 


•  i 


G.  . 

•J 


P' 

'j 


(4-52) 
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for  the  case  where 
equation. 


f  >  0,  the  stress  rate  may  be  expressed  by  the  following 


If  time  is  considered  the 
the  incremental  stress  is 


integration  parameter  and  the  time  step  is 


At, 


do 


ij 


AdGkk6i j 


The  absence  of  plastic  volume  strain 
rule  in  which  Jj  does  not  appear  In 


s  due  to  choosing  a  nonassoc iated  flow 
the  plastic  potential . 


S°me  example  calculations,  which  are  summarized  in  Tables  4-3  and  4.4 
are  shown  in  Fibres  4-11  throu9h  4-13.  Comparison  with  some  data  for  a 
tonal, te  ,s  shown  in  Figure  4-14.  The  comparison  illustrates  the  ability  of 
t  e  model  to  represent  increase  ir  strength  with  strain  rate  and  its  inability 
to  represent  nonlinear  behavior  prior  to  yielding.  1,  can  be  shown  that  the 
e  feet,  of  v.scoplasticity  can  be  accounted  for  entirely  by  a  correction  ,0  the 

I °ad  W'Ct°r  (Refer'"Ce  "-W-  Th'  technique  is  similar  to  that  described  above 
for  rate- Independent  plasticity  and  to  that  described  below  for  viscoelasticity. 

An  important  consequence  of  this  is  that  time-consuming  reformulation  of  the 
stiffness  matrix  is  avoided. 
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TABLE  4-3.  SUMMARY  OF  EXAMPLE  CALCULATIONS 


Case 

Type  of  Loading 

Stress  Ratio, 
°3/ol 

Strain  Rate,* 
in. /in. -sec 

1 

—  —  — 

Proportional  loading 

0.293 

200.00 

2 

Proportional  loading 

0.293 

2.00 

3 

Proportional  loading 

0.293 

0.20 

4 

Proportional  loading 

0.293 

0.0? 

5 

Proportional  loading 

0.200 

0.20 

6 

Proportional  loading 

0.100 

0.20 

7 

Uniaxial  strain 

0.293 

0.02 

8 

Uniaxial  strain 

0.293 

0.20 

9 

Uniaxial  strain 

0.293 

2.00 

10 

Uniaxial  strain 

0.293 

200.00 

’'For  proportional  loading,  these  are  elastic  strain  rates. 


TABLE  .  PROPERTIES  USED  IN  PRESENT  EXAMPLES 
Propert ies 


Bui k  modulus 

B 

=  1.5x10 

Shear  modulus 

G 

=  10^  psi 

Cohesion 

a 

0 

=  1450 

Angle  of  friction 

al 

=  -0.1 

Viscoplastic  coefficient 

Y 

=  1  .0 

FIGURE  4-13.  COMPRESSIVE  STRESS/STRAIN  RELATIONS 
FOR  PRESENT  MODEL  SUBJECTED  TO 
PROPORTIONAL  LOADING 


FIGURE  4-14.  COMPRESSIVE  STRESS/STRAIN  RELATIONS  FOR  CEDAR  CITY  TONAL  I TE  AT 
VARIOUS  STRAIN  RATES 
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*♦.1.6  RATE  EFFECTS— VISCOELASTIC iTY 

The  total  strain  is  defined  to  be  the  sum  of  instantaneous  elastic 
and  viscoelastic  parts.  The  strain  is  further  divided  into  shear  and  volumetric 
components,  which  are  treated  separately  as  follows: 


£«•  —  5..  e,  . 

ij  ij  kk 


(4-55a) 


+  (ekk*C 


(4-55b) 


In  the  computer  program  the  user  may  choose  to  have  either  elastic  or  visco¬ 
elastic  shear  deformation  and  to  have  either  elastic  or  viscoelastic  volumetric 
deformation. 

Kelvin,  Maxwell  and  three-parameter  fluid  models  are  available  as 
shown  in  Figure  4-15.  To  simplify  the  following  discussion,  no  distinction  is 
made  between  volumetric  and  shear  components.  The  creep  rate  and  creep  strain 
at  t  +  At  may  be  expressed  as  follows: 


KELVIN 


•c  .  c 
e  +  a1  e  =  a2o 


(4-56a) 


+  At  =  et  6Xp('alAt)  +at  ^(1  '  eXp('al' 


(4-56b) 


E  ,  1 

where  a^  =  —  and  a2  =  —  ,  E  and  n  are  the  spring  and  dashpot  constants 
respectively. 


MAXWELL 


e  -  a^  +  a2a 


st  +  At  =  Et  +  [(a2At  +  a1}  °t  + 


(4-57a) 


(4-57b) 


where  a,  =  i  and  a_  =  — 
IE  2  n 


A 
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KELVIN 


3 


“WV 


MAXWELL 


3 


3  PARAMETER 
FLUID 


FIGURE  4-15.  VISCOUS  MODELS  AVAILABLE  IN  THE 
PRESENT  COMPUTER  PROGRAM 
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three-parameter  fluid 


The  creep  strain  of  the  three-parameter  fluid  is  divided  into  two 

components  e  and  e  which  represent  the  creep  strain  in  a  Kelvin  model 
and  dashpot 


c  _  1  c  .  2  c 

e  -  e  +  e 


(4-58a) 


1  c 

Gt  +  At 


'et  exP(_a1At)  +  ot  —  (1  "  expf-a,  At)) 


(4-58b) 


2  c 

Gt  +  At 


2  c  4.  At  * 

£t  +  At  +  57  ot 


(4-58c) 


These  models  are  illustrated  in  Figure  4-15.  Various  aspects  of  their  proper¬ 
ties  are  discussed  in  Reference  4-35.  Equations  4-56a,  -57a,  and  -58a  are  in 
suitable  form  for  application  to  a  time-marching  Integration  procedure. 

Typical  creep  and  creep  recovery  data  for  several  rocks  are  shown  in 
Figure  4-16.  Application  of  a  Kelvin  model  to  one  of  the  rocks  is  shown  in 
Figure  4-17.  By  suitably  varying  the  parameters  a?  and  a  ,  a  range  of 
behavior  can  be  reproduced.  To  assist  the  inexperienced  user  in  selective 
coefficients  for  these  models,  the  following  guidelines  are  offered. 

The  finite  element  adaptation  of  these  models  uses  the  initial 
strain  approach  to  writing  the  equilibrium  equations  based  on  the  change  in 
internal  energy. 


JaTe  dV  -  f  PTu 
vo  1  J  S 


(4-59) 


The  stress/strain  relations  are  in  terms  of  a  matrix  C  of  elastic  stress/ 
strain  coefficients 


V  -  lC] }eQ  +  dc  -  e£j 


(4-60) 
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(a)  CREEP  AND  CREEP  RECOVERY  (SLATE)— REFERENCE  4-28 


(b)  CREEP  UNDER  CONSTANT  STRESS  AND  RECOVERY  CURVES --REFERENCE  4-30 
(AFTER  REFERENCE  4-31) 


FIGURE  4-16.  EXPERIMENTAL  DATA  ON  VISCOELASTIC  PROPERTIES  OF  SEVERAL  ROCKS 
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COMPARISON  BETWEEN  PRESENT  CREEP  MODEL  AND  EXPERIMENTAL  DATA 
(REFERENCES  J»-30,  J»-31) 
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where 


Stress  at  t 

Total  strain  (including  creep)  at  t  =  t 

Strain  increment  (including  creep)  during  At  =  t  -  tQ 

Accumulated  creep  to  t  =  t 


An  error  is  committed  in  Equation  4-60  in  that  de  contains  both  elastic  and 
viscous  components,  whereas  it  is  treated  as  if  it  were  entirely  elastic. 
Defining  the  strain/displacement  transformation  matrix  B  as 


{e}  =  [ B] {u } 


(4-61) 


The  strain  energy  for  an  element  is 


U  = 


BTCBdV  u  + 


(4-62) 


Def i n i ng 


BTCBdV  =  k 


(4-63) 


where  K  is  the  element  stiffness  matrix  and  performing  a  variation  with 
respect  to  the  generalized  displacements  u  results  in 

kdu  =  P  -  F  (4-64) 

e 

where 
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It  is  significant  that  the  element  stiffness’ matrix  consists  entirely  of 
linear  elastic  terms.  Thus,  the  assemblage  stiffness  matrix  needs  to  be 
formulated  only  once,  resulting  in  great  economy  of  computing.  It  is  neces¬ 
sary  to  store  the  components  of  the  creep  strain  at  t  for  use  in  the  next 

o 

step. 

This  approach  has  been  adapted  to  finite  element  for  rock  and 
concrete  (Reference  -1-26) .  Some  experimental  work  which  has  been  performed 
on  rock  is  compared  in  Reference  4-27  with  a  series  of  Kelvin  models. 

Reference  4-28  discusses  methods  of  accounting  separately  for  volumetric  and 
deviatoric  creep  strains. 

4.2  MATERIAL  PROPERTIES  OF  JOINTS 

This  section  describes  the  properties  which  may  be  assigned  to  joints. 
These  properties  consist  of  the  shearing  and  normal  stiffnesses  of  the  joints. 
They  correspond  physically  to  the  stiffness  and  strength  of  fault  gouge,  to 
the  roughness  of  the  joints  and  to  the  angles  of  slip  surfaces  relative  to  the 
principal  plane  of  the  joint.  They  are  classed  as  dilatant  if  shearing  pro¬ 
duces  joint  expansion  or  contraction  or  nondilatant  if  shearing  and  normal 
displacement  are  uncoupled.  The  properties  are  specified  by  the  user  in 
natural  coordinates  which  may  be  directions  parallel  and  perpendicular  either 
to  slip  surfaces  or  to  the  principal  plane  of  the  joint.  In  either  case, 
transformation  to  global  directions  is  automatically  performed  by  the  program. 

As  in  the  case  of  homogeneous  material  properties,  the  joint  prop¬ 
erties  are  controlled  by  subroutines  CONECT  and  ELPL  through  modular  subroutines. 
Presently,  these  contain  built-in  joint  properties.  As  more  data  become 
available  on  properties  of  joints  and  the  present  joint  material  properties 
become  obsolete,  the  present  model  may  easily  be  modified  without  disturbing 
the  main  program  or  any  of  thf  other  material  properties. 
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4.2.1  NONDILATANT  JOINTS 

This  class  of  Joints  is  the  simplest  to  model  mathematically  since 
ere  ,s  no  volume  change  due  to  shearing  strains,  and  therefore  the  shear 

and  the  normal  components  of  deformation  are  uncoupled  and  the  stress-strain 
relations  are  as  follows: 


C  0 
nn 


(4-65) 


However,  Cnn  and  Cgs  are  nonlinear  functions. 

In  stress-deformation  relationship  in  normal  direction,  three 
istinct  stages  can  be  recognized  (see  Figure  4-18); 

a.  Separation,  Cnn  =  =  0  when  ^  >  0 

b.  Crushing  of  the  surface  irregularities  or  the  compression  of 
the  material  in  the  fault  or  joint,  if  any  C  =  E 

(en  «  en  "  0)’  For  smooth  surfaces  this  case"does  not  exist, 
therefore  eC  =  0. 


c.  Contact,  C  =  E  (e  <  ec) 
nn  r  n  n 


It  Is  important  to  note  that  very  high  values  can  be  assigned  to  E,  without 
any  numerical  problems  with  the  special  join,  element  described  Section  3. 

The  tangential  stress-strain  relationship  Is  assumed  to  be  elastic- 
perfectly  plastic  using  a  Mohr-Coulomb  yield  criterion: 


C  =  G 
ss 


C  =0 
ss 


°c  <  c  +  o  tan 
s  n 


o  -  c  +  o  tan 
s  n 


(4-66) 


where  c  and  <j>  are  the  cohesion  or  the  angle  of  frictic 
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FIGURE  4-18.  NORMAL  STRESS-STRAIN  RELATION  FOR  JOINTS 
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4.2.2  DILATANT  JOINTS 

Bilatancy  of  rock  joints  and  faults  are  very  complex  to  model 
mathematically;  however,  to  include  a  measure  of  dilatancy,  the  procedure 
developed  by  Goodman,  Dubois,  and  Brekke  (Reference  4-36)  is  used  here. 

Further  data  are  available  from  Reference  4-37. 


It  is  assumed  that  the  deformation  in  p-q  coordinate  system  shown 
in  Figure  4-19  is  nondilatant.  The  angle  y  between  the  p  direction  and  the 
j~int  surface  is,  therefore,  defined  to  be  a  material  property  of  the  joint. 

The  stress-strain  relation  in  n-s  coordinate  is: 


4.2.3  DEMONSTRATION  OF  JOINT  ELEMENT 

The  example  problem  in  Figure  4-20  demonstrates  the  behavior  of  the 
present  joint  element.  The  problem  consists  of  a  plane  system  of  three  solid 
blocks  and  three  joint  planes.  The  system  is  restrained  at  the  bottom  and  is 
subjected  to  lateral  pressure  P  and  vertical  pressure  2P.  The  joints  are 
assumed  to  be  filled  with  nondilatent  material  having  the  properties  shown  in 
Figure  4-20,  The  blocks  are  assumed  to  have  the  isotropic  elastic  properties 
shown  in  Figure  4-20.  The  material  in  the  joints  is  represented  by  joint 
elements  which  are  shown  as  their  shaded  strips.  Ambiguity  with  respect  to 
nodal  displacements  would  arise  at  the  point  of  the  wedge-shaped  block  if  the 
joint  elements  were  to  extend  to  the  point.  Ambiguity  is  avoided  by  stopping 
the  elements  short  of  the  point.  This  is  a  good  approximation  to  a  physical 


I 

! 
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geometry  of  example  problem 
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system  if  the  joint  elements  extend  infinitesimally  close  to  the  point.  In  the 
present  example,  the  gap  is  made  large  to  Illustrate  how  the  difficulty  is 
avoided.  The  blocks  are  represented  by  plane  strain  elements. 

Results  of  the  present  example  are  shown  in  Figures  4-21  and  4-22. 
Figure  4-21  shows  how  the  vertical  displacement  varies  as  a  function  of  hori¬ 
zontal  distance  from  the  centerline.  Each  curve  in  Figure  4-21  corresponds 
to  a  different  distance  below  the  top  surface.  Figure  4-22  shows  the  wedge  in 
its  final  position. 


A 


FIGURE 
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-22.  FINITE  ELEMENT  MESH  IN  THE  DISPLACED  CONFIGURATION  FOR 
WEDGE  PROBLEM  IN  FINAL  EQUILIBRIUM  STATE 
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SECTION  5 

DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

This  section  describes  the  computer  program  in  terms  of  features 
which  are  apparent  to  the  user,  such  as  automatic  mesh  generation,  and  in 
terms  of  its  logical  structure.  The  first  part  of  the  section  describes  the 
steps  which  the  user  takes  in  order  to  prepare  the  input  data.  The  mesh 
generator,  bandwidth  reducer  and  plotting  capability  are  described.  The 
motive  and  logic  diagrams  for  element  renumbering  are  also  described. 

The  second  part  of  this  section  describes  the  execution  part  of  the 
computer  program  in  terms  of  logic  diagrams.  Subroutines  are  described  which 
control  the  computation,  form  the  global  stiffness  matrix  and  compute  the  load 
vector,  form  element  stiffnesses  and  other  operations.  Also  described  is  the 
technique  of  multibuffering  whereby  data  is  transferred  between  core  and  per¬ 
ipheral  storage  units  at  the  same  time  that  computations  are  being  performed 
in  core.  This  technique  greatly  improves  the  efficiency  of  the  program  for 
problems  where  large  blocks  of  data  are  stored  on  peripheral  memory  units. 

5.1  PROCESSING  OF  INPUT  DATA 

The  sequence  in  which  input  data  are  processed  is  shown  in  Figure  5_1 • 
There  are  three  basic  phases.  The  first  consists  of  generating  continuum 
(two-  or  three-dimensional)  and  joint  elements  and  plotting  the  results.  The 
second  phase  consists  of  adding  other  elements  (beam,  thick  shell,  bar,  etc.) 
manually,  reducing  the  bandwidth,  shuffling  the  element  numbers  so  that  they 
appear  in  the  order  of  contribution  to  the  global  stiffness  matrix,  and 
plotting  the  final  mesh.  The  third  phase  consists  of  reading  additional  data 
such  as  material  properties  and  loads. 

MESH  GENERATOR 

The  automatic  mesh  generation  scheme  incorporated  in  the  program 
requires  the  user  to  provide  a  coarse  mesh  which  the  program  refines  by 
subdivision  under  the  control  of  the  user.  The  main  goal  is  to  minimize  the 
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FIGURE  5-1.  OPERATIONS  OF  THE  INPUT  PORTION  OF  THE  PRESENT  COMPUTER  PROGRAM 
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the  input  preparation  on  the  part  of  the  user.  The  method  uses  the  key 
diagram  concept  described  in  Reference  5-1;  the  remainder  of  the  present 
mesh  generation  scheme  is  new  and  original  work. 

The  subdivision  of  the  given,  coarse  mesh  to  obtain  the  final,  fine 
mesh  is  carried  out  by  subdividing  each  one  of  the  mesh  units,  as  "zones," 
within  the  coarse  mesh  in  the  following  manner.  Consider  the  simple  case 
or  a  general,  two-dimensional  mesh,  in  which  the  basic  mesh  unit  or  zone  is 
a  parabolic  quadrilateral  as  shown  in  Figure  5-2.  Assume  that  the  x  and  y 
coordinates  of  the  four  corner  mesh  points,  1,  2,  3,  and  4,  as  well  as  those 
of  the  four  m>dpoints,  5,  6,  7,  and  8,  are  known.  The  x  and  y  coordinates 
of  an  arbitrary  point  P  within  this  zone  can  then  be  expressed  conveniently 
in  terms  of  the  coordinates  of  the  eight  points  through  the  use  of  the  local, 
curvilinear  coordinates,  £  and  n>  whose  values  range  from  -1  to  +1  on 
opposite  sides  as  indicated  in  Figure  5-2.  Thus,  one  can  write 


8 

x  =  y.  n.x. 

M  1  1 

8  (5" 

y  =  E  N.y. 

i=i 

in  which  N.'s  are  the  so-called  "isoparametric  shape  functions"  expressible 
in  terms  of  the  curvilinear  coordinates  as  follows: 

N!  "  -  0  -  OO  -  n)U  +  n  +  D 

N2  "  "  ^  ^  "  0(1  +  n)  U  -  n  +  1) 

N3  "  T  0  +  0  (1  +  n)  U  +  n  -  1) 

n4  ■  £  (i  +  0(1  -  n)U  -  n  -  D 

N,  -  j  (1  -  5)0  -  n2)  (5’21 

n6  -  \  (i  -  Od  -  n2) 

N?  -  7  0  -  e2)0  +  n) 

N8  •  7  0  -  S2)(i  +  n) 
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FIGURE  5-2. 


QUADRILATERAL  ZONE,  WHOSE  SHAPE  AND  COORDINATES  ARE 
EXPRESSED  BY  PARABOLIC  SHAPE  FUNCTIONS 
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Thus,  for  a  certain  set  of  values  of  the  curvilinear  coordinates,  or  "natural 
coordinates,"  the  corresponding  Cartesian  coordinates  can  be  easily  found 
from  Equations  5-2.  Consequently,  to  subdivide  the  zone  into  a  A  x  A  mesh, 
for  example,  and  establish  the  Cartesian  coordinates  of  the  twenty-one  new 
mesh  points,  one  merely  substitutes  into  Equations  5"1  twenty-one  times, 
each  time  with  a  different  combination  of  values  of  e  and  n  which  are 
incremented  successively  by  0.5. 

This  technique  can  be  used  to  treat  each  one  of  the  zones  in  a 
given,  coarse  mesh  of  any  configuration.  The  only  restriction  that  has  to  be 
observed  is  that  the  number  of  subdivisions  between  any  two  adjoining  zones 
must  match.  This  is  to  satisfy  the  fundamental  connectivity  requ i rt.T^nt  of 
the  finite  element  method.  To  satisfy  this  continuity  requirement  and  to 
facilitate  the  preparation  of  the  input  that  defines  the  basic  coarse  mesh, 
it  is  convenient  to  introduce  the  so-called  "key  diagram." 

A  key  diagram,  in  general,  is  a  rectangular  grid  resembling  a 
checker-board.  It  has  no  physical  dimensions.  Its  purpose  is  to  present, 
or  define,  the  connectivity  of  the  zones  in  the  coarse  mesh  and  to  facilitate 
defining  the  position  of  a  mesh  point  in  the  form  of  row  and  column  numbers 
in  the  coarse  mesh.  Another  purpose  of  the  key  diagram  is  to  help  define  the 
extent  of  subdivision  of  a  zone,  for,  to  satisfy  the  aforementioned  connec¬ 
tivity,  a  newly  introduced  mesh  line  has  to  extend  across  all  the  zones 
located  in  the  same  row  or  column  in  the  key  diagram. 

Figure  5“3  illustrates  the  generation  of  a  simple  finite  element 
mesh  representing  a  dam  and  part  of  its  foundation.  The  domain  has  been 
blocked  into  three  zones  which  are  connected  as  shown  in  the  accompanying 
key  diagram.  The  final  mesh  is  the  result  of  subdividing  the  first  and  the 
last  zone  into  two  subdivisions  and  the  second  into  three  subdivisions  in 
the  vertical  direction,  and  all  three  zones  into  three  subdivisions  in  the 
horizontal  direction. 
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of  .he  M2,T'  'he  Pr°9ram  f°r  th'S  PartiC“'ar  eXa"P,e  *U,d 

a-  A  statement  .ha.  the  Key  diagram,  or  the  origins,  coarse  mesh, 
• S  3  by  1 .  * 


b. 


The  rectangular  coordinates  of  the  given,  eight  mesh  points 
and  those  of  the  single  midpoint  that  define  the  lower  boundary 

OH  A  mrwA  * 


on  a  curve. 

c.  Three  materia,  numbers  to  associate  the  firs,  two  zones  with 
concrete  and  the  third  with  the  soil  medium. 

"■  rfThe.‘hrae  3  ^  2>  ‘hat  SP<!ciff  the  number  of  sub- 

divisions  in  the  vertical  direction,  and  the  number  2  for 
the  horizontal  direction. 

The  computer  program  would  assign  the  materia,  number  for  a  zone  to  all  of 
the  elements  that  are  created  In  -* 

elements  in  rh.  f  ,  !  ’  "d  the  noda'  *>'»«»  and 

left  ,  j  '  81  meSh  f™»  th.  top  to  the  bottom,  from  the 

, toward  the  right  by  referring  to  the  Key  diagram.  The  curved  mesh  Has 

*  t  H  “ne  are  meant  *»  the  fact  that  nodal  points  are 

senerated  on  a  second  degree  curve  because  of  the  nature  of  the  shape  func 

r^T'  ma!h  . . -  a  -ight  ,,„a  between  two!:::;I:;ng 

In  the  previous  example.  If  the  user  Is  not  satisfied  with  the 
re  inement  of  the  final  mesh  or  with  that  in  any  one  of  the  .  ca„ 

qu.cKly  ,eru„  the  problem  by  modifying  the  four  Integers  mentioned  „em  d 

w i thou t  having  to  change  any  of  the  data  mentioned  In  the  first  three  I. 

Wh;;;  .fin  he  given,  coarse  mesh.  This  Is  perhaps  the 
re  o  this  particular  mesh  generation  scheme. 
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Aside  from  being  used  to  define  an  edge  of  a  zone  on  a  curve,  a 
midpoint  can  also  be  used  to  imply  a  mesh  grading.  The  mesh  shown  in 
Figure  is  the  result  of  specifying  the  three  midpoints,  M^ ,  M2  and  M^, 
purposely  off-centered  and  toward  the  edge  of  the  circular  hole.  Midpoints 
M^  and  M2  are  there  merely  to  form  the  edge  of  the  hole  on  a  curve.  In 
this  particular  example,  number  of  subdivisions  of  5  and  5  in  the  vertical 
direction  (in  the  key  diagram)  and  of  16  in  the  horizontal  direction  were 
speci f ied. 


Another  feature  that  makes  the  scheme  versatile  and  powerful  is 
the  provision  for  the  user  to  prescribe  a  zone  as  void  merely  by  assigning 
zero  for  the  material  number  of  that  zone.  With  this  provision,  one  can 
easily  generate  a  mesh  around  a  notch  or  a  cutout.  The  mesh  around  a  tunnel 
shown  in  Figure  5"5  is  such  an  example. 

The  final  feature  to  be  mentioned  here  is  the  program's  ability  to 
join  any  two  edges  in  the  key  diagram.  This  makes  it  possible  to  produce  the 
meshes  shown  in  Figures  5_6  and  5“7<  The  user  may  prefer  the  type  of  mesh 
shown  in  Figure  5"7  over  that  shown  in  Figure  5"5  for  the  radiating  pattern 
of  the  mesh  lines  which  makes  it  easier  to  make  a  finer  mesh  around  the 
edge  of  the  tunnel . 

The  same  scheme  can  be  used  to  generate  a  three-dimensional  mesh. 
The  necessary  changes  are  to  replace  the  two-dimensional  shape  function  and 
key  diagram  by  their  three-dimensional  counterparts,  and  to  introduce  the 
third  Cartesian  coordinate.  The  basic  "zone"  or  "block"  in  this  case  is  a 
general,  six-sided  body  with  parabolic  edges.  It  has  eight  corner  nodes  and 
twelve  midpoints  as  shown  in  Figure  5~ 8 ,  and  the  Cartesian  coordinates  at  an 
arbitrary  point  P  within  the  block  having  the  natural  coordinates  £,  n 
and  z  are  expressible  in  terms  of  the  Cartesian  coordinates  of  the  twenty 
controlling  points  as  follows: 
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(5-3) 


where,  for  corner  points, 

N.  »  |  (l  +  CjC)(1  +  n;n)  (1  +  C .  c)  ^  +  n.n  +  £.£  -  2)  (5-M 

and  for  midpoints, 

N.  -  ^  (1  -  £2)(1  +  n.n)(l  +  £jO  for  points  with  £{  ■  0 

N.  -  ^  (1  +  S.OO  -  n2)(l  ♦  CjC)  for  points  with  r1{  -  0 

N .  -  ^  (1  +  £.0(1  +  n.n)(l  -  £2)  for  paints  with  £.  -  0 

As  before,  zero  can  be  used  as  material  number  to  imply  void  or  cavity,  and 
a  midpoint  can  be  so  specified  as  to  imply  a  curved  edge  as  well  as  a  mesh 
grading  along  the  edge.  Instead  of  joining  two  edges,  one  can  now/  join  any 
two  interfaces  as  long  as  the  meshes  on  them  are  compatible.  Some  three- 
dimensional  meshes  automatically  generated  in  this  manner  are  illustrated 
in  Figures  5-9  and  5*10. 
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BANDWIDTH  REDUCER 

The  computer  time  tequlred  to  perform  a  computation  is  approximately 
proportional  to  the  square  of  the  bandwidth.  Thus  It  is  Important  to  aim 
for  a  minimum  bandwidth  when  numbering  nodal  points.  However,  when  the  finite 
element  mesh  is  generated  automatically  py  the  technique  described  above, 
no  attention  whatever  is  paid  to  minimizing  bandwidth.  The  user  can  influence 
the  bandwidth  to  some  extent  by  judicious  choice  of  key  diagram.  Neverthe¬ 
less,  it  Is  very  iikely  that  the  bandwidth  so  generated  will  not  be  optimum. 
The  situat  on  wi i i  often  be  made  worse  when  elements  are  added  manually. 

To  help  alleviate  this  difficulty  and  thus  encourage  the  user  to  use  the 
automatic  mesh  generator,  a  bandwidth  reducer  is  Included  in  the  present 
computer  program. 

The  function  of  the  bandwidth  reducer  is  shown  In  Figure  5-11.  The 
INPUT  NODAL  POINT  configuration  is  typical  of  the  mesh  which  would  be  gen¬ 
erated  automatically  for  which  the  maximum  difference  in  nodal  point  numbers 
for  any  element  Is  31.  This  configuration  was  submitted  to  the  bandwidth 
reducer.  The  result,  after  2  sec  of  computation  time  by  the  bandwidth 
reducer  is  the  REDUCED  NODAL  POINT  configuration,  for  which  the  maximum 
difference  In  node  numbers  Is  6.  That  this  configuration  is  not  optimum  is 
shown  by  the  IDEAL  configuration  for  which  the  maximum  difference  is  5.  Thus, 
the  technique  is  not  optimum  because  it  does  not  always  converge  to  the 
minimum  bandwidth  in  the  computer  time  which  the  user  has  selected.  Aiso, 
the  technique  operates  on  nodal  point  numbers  rather  than  on  degree  of 
freedom  numbers.  Thus,  when  several  types  of  elements  are  mixed,  the  con¬ 
figuration  corresponding  to  minimum  difference  in  node  numbers  does  not 
necessarily  correspond  to  minimum  bandwidth. 

Logic  diagrams  for  the  bandwidth  reducer  are  shown  in  Figures  5-12 

and  5-13- 

ELEMENT  NUMBERING 

In  addition  to  renumbering  nodal  points  to  reduce  the  bandwidth  of 
the  global  stiffness  matrix,  the  elements  are  numbered  in  the  order  cf  their 
contribution  to  the  global  stiffness  matrix.  In  this  way,  the  strain/ 
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INPUT  NODAL  POINT  CONFIGURATION 
MAXIMUM  NODAL  POINT  DIFFERENCE  -  3< 


REDUCED  NODAL  POINT  CONFIGURATION 
MAXIMUM  NODAL  POINT  DIFFERENCE  -  6 
SOLUTION  TIME  -  2  SECONDS 


IDEAL  CONFIGURATION 

MAXIMUM  NODAL  POINT  DIFFERENCE  -  5 


AAM8I 


FIGURE  5- II.  EXAMPLE  OF  BANDWIDTH  REDUCER 
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DETERMINE  THE  MAXIMUM  BANDWIDTH  AND  THE  TWO 
NODES  CAUSING  THE  MAXIMUM  BANDWIDTH 


CAN  THE  HIGHER  NUMBERED^***. 
NODE  BE  INTERCHANGED  WITH  A 
LOWERED  NUMBERED  ONE  TO  REDUCE 
THE  BANDWI DTH? 


CAN  THE  LOWERED  NUMBERED*'** 
NODE  BE  INTERCHANGED  WITH  A 
HIGHER  NUMBERED  ONE  TO  RE¬ 
DUCE  THE  BANDWIDTH?  ^ 


""CAN  THE  HIGHER  NUMBERED  NOD?* 
BE  INTERCHANGED  WITH  A  LOWER 
NUMBERED  NODE  AND  STILL 
.MAINTAIN  THE  SAME  BANUWIDTHL 


CAN  THE  LOWER  NUMBERED  NOdTn>>^ 
BE  INTERCHANGED  WITH  A  HIGHER 
NUMBERED  NODE  AND  STILL  MAINTAIN, 
.JHE  SAME  BANDWI  OTH? 


NO  FURTHER  INTERCHANGES  ARE  POSSIBLE 
EXIT 


PERFORM  THE  INDICATED  NODE  INTERCHANGE 


MODIFY  TME  "IS"  TABLE  TO  REFLECT  THE  CURRENT 
INTERCHANGE 


AAA6S9 


FIGURE  5-13.  PROCEDURE  USED  IN  BANDWIDTH  REDUCER  (REFERENCE  5-3) 
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A 

stiffness  matrixes  can  be  retrieved  from  peripheral  storage  in  the  least 
time.  Since  the  element  data  are  used  in  the  basic  operations  of  forming 
the  giobal  stiffness  matrix  and  effective  load  vector,  the  efficiency  gained 
oy  performing  these  complicated  operations  Is  essential. 

The  element  numbers  assigned  by  the  mesh  generator  and  the  user 
are  revised  such  that 

min(LM(l ) )  for  the  nth  element  < 
m i n ( LM ( I ) )  for  the  n  +  1^  element 


where 

LM(l)  -  Array  of  degree  of  freedom  members  for  the  element 

Element  data  are  stored  sequentially  on  peripheral  storage  units  such  that 
the  data  for  Element  1  is  at  the  head  of  the  unit,  and  then  in  sequence 

Element  Number  -  1  <  ...n,n+1...  <  NUMEL 

where 

NUMEL  ■  Total  number  of  elements 

In  order  to  modify  the  effective  load  vector  F,  to  obtain  the 
matrix  C  of  stress/strain  coefficients  and  to  obtain  element  stresses  and 
strains,  an  array  containing  the  previous  loads  F,  the  displacements  u  and 
the  incremental  displacements  du  is  brought  into  core  as  shown  in  Figure  5-15. 
Also  required  are  data  stored  on  the  element  data  tape  such  as  the  strain/ 
displacement  transformation  matrixes  and  LM  arrays  for  the  elements  of 
interest.  As  shown  in  Figure  5-14,  all  degrees  of  freedom  to  which  Element  n 
contributes  are  found  in  a  sequence  of  F,  u,  du  array  one  bandwidth  (MBAND) 
long.  Thus  it  is  possible  to  process  to  F,  u,  du  jirray  by  reading  the 
element  data  tape  only  once. 
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is  added  to  Block  nn.  Now  Block  nn  is  completely  prepared  and  is  moved  out 
of  core  while  Block  nn  +  1  is  brought  into  core.  Block  nn  +  1  is  filled  in 
the  same  way  as  Blocks  nn  -  1  and  nn,  which  is  by  passing  sequentially  through 
the  element  data  array.  When  all  direct  contributions  to  Block  nn  +  1  have 
been  made  and  all  data  which  contributes  to  Block  nn  +  2  has  been  stored  on 
unit  T2,  the  data  stored  on  unit  T1  is  read  into  core  and  added  to 
Block  nn  +  1 .  This  block  is  now  completely  prepared  and  is  transferred  to 
peripheral  storage. 

SIMULATION  OF  THE  CONSTRUCTION  AND  EXCAVATION  SEQUENCE 

The  construction  and  excavation  sequence  will  be  represented  by 
modification  of  the  material  properties  of  the  elements  involved.  Initially 
elements  will  be  assigned  to  the  regions  to  be  constructed  or  excavated  later. 

A  flag  will  indicate  the  time  of  the  addition  or  subtraction  of  each  element. 

In  case  of  excavation,  the  elements  of  that  region  will  have  the  appropriate 
material  properties  and  will  contribute  to  the  global  stiffness  of  the  system 
up  to  the  indicated  time,  beyond  which  the  contribution  of  these  elements  to 
the  global  stiffness  of  the  system  will  be  zero.  The  reverse  of  this  procedure 
will  be  applied  to  the  elements  of  the  regions  to  be  constructed  later.  All 
these  operations  will  be  performed  in  the  element  package. 

SELF-LOADING 

The  initial  state  of  the  system  will  be  assumed  to  be  stress  free. 
Then  the  dead  load,  if  any,  will  be  applied  in  a  specified  number  of  increments 
prior  to  the  application  of  the  external  load.  It  is  necessary  to  apply  the 

self-loading  incrementally  due  to  the  fact  that  the  system  is,  in  general, 
nonl i near. 
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5.2  EXECUTION  PHASE  OF  PROGPAM 

The  execution  phase  of  the  program  is  summarized  in  Figures  5-16 
through  5-21.  The  controlling  program  is  BMCALC ,  Figure  5-1 6 .  The  first 
main  operation  is  to  form  the  effective  load  vector  and  global  stiffness 
matrix,  which  is  performed  by  Subroutine  KFORM ,  Figure  5-17.  The  assembly 
of  the  global  stiffness  from  the  element  stiffnesses  is  performed  by  Sub¬ 
routine  BSTIF,  Figure  5-18.  Subroutines  FWRT ,  TDRUM  and  FILLFU  transfer  data 
from  peripheral  storage  to  core. 


5.2.1  MULTIBUFFERING  TECHNIQUE 

Multibuffering  is  a  technique  whereby  central  processor  wait  time 
for  all  binary  Read/Write  operations  involving  the  peripheral  storape  of  data 
is  minimized.  (Formatted  I/O  operations  such  as  card  reading,  punching,  and 
printing  are  not  included  in  this  discussion.)  The  reason  for  using  multibuf¬ 
fering  techniques  is  that  data  moves  faster  between  core  locations  than  between 
peripheral  storage  and  core. 


The  problem  which  mu  1 1 i buf fer ing  overcomes  is  the  standard  I/O 
feature  of  higher  level  programming  languages,  such  as  FORTRAN,  which  requires 
that  when  an  I/O  operation  (READ  or  WRITE)  is  initiated,  computation  ceases 
until  the  I/O  operation  is  completed.  This  feature  assures  the  user  that  data 
he  may  wish  to  use  is  in  core  before  he  tries  to  use  it.  The  amount  of  time 
the  program  must  wait  for  completion  of  I/O  operations  depends  on  (1)  the 
access  time  and  (2)  the  data  transfer  time.  These  times  depend  on  the  type 
of  peripheral  device  being  used  and  the  amount  of  data  to  be  transferred. 

Thus,  as  the  amount  of  peripheral  storage  increases,  the  time  spent  waiting 
for  completion  of  I/O  operations  increases  and  for  large  volumes  of  data 
the  I/O  time  may  control  overall  run  times.  Mu  1 1 i buf fer i ng  minimizes  the 
wait  time  of  these  I/O  operations  by  allowing  computations  to  proceed  at  the 
same  time  data  transfer  from  peripheral  storage  is  occurring.  This  requires 
standard  FORTRAN  I/O  operations  to  be  replaced.  This  is  possible  on  most 
large  scale  scientific  computers  by  the  use  of  special  machine-dependent 
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FIGURE  5-16.  PROGRAM  BMCALC— CONTROLS  MAIN  OPERATIONS  OF  THE 
COMPUTATION  SECTION 
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CALL  FWRT  TO  HOVE 
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FIGURt  5 “17.  (CONTINUED) 
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FIGURE  5-19. 


SUBROUTINE  FWRT— MOVES  DATA  FROM  LIVE  LOAD  VECTOR  (FVORK) 
TO  F,  u,  du  uiiTPUT  BUFFER  AND  TO  u  OUTPUT  FILE 
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FIGURE  5-20. 
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statements  or  -.routines.  These  ,„terfece  routines  ere  generally  different 
on  oech  machine  end,  unless  cere  Is  exercised,  e  program  may  become  machine 
ependent.  This  danger  Is  avoided  In  the  present  application  of  multlbuff.r- 
"9  y  Isolating  all  I/O  statements  In  one  subroutine  which  may  easily  be 
modified  when  moving  the  program  to  a  different  machine. 

A  general  multibuffering  scheme  to  perform  computations  on  a  set  of 
data  stored  on  peripheral  devices  Is  shown  In  Figure  5-22.  This  scheme, 

Which  Is  Incorporated  In  the  present  computer  program,  requires  either  the 
amount  of  main  storage  used  for  buffers  to  be  Increased  or  the  size  of  each 

b'  d'Cr'asad  relative  to  the  buffer  size  that  might  be  used  with 
standard  FORTRAN  I/O  procedures.  Since  most  of  the  main  storage  Is  already 

test  In"  d'"rm'n'n9  bUff8r  tha  ,at"r  a'«rna,  Ive  Is  employed.  Extensive 

testing  and  previous  experience  have  shown  that,  although  the  number  of  I/O 

operations  Increases,  multibuffering  results  In  substantial  overall  reduction 

n  onputer  run  time  for  a  given  problem.  Savings  Increase  with  the  problem 

size  and  thus  the  volume  of  data  on  peripheral  storage  Increases. 

5.2.2  BAND  SOLVER 

The  solution  of  large  structural  systems  requires  the  solution  of 
a  set  of  linear  simultaneous  equations  of  the  form 


-  [K]  {u} 

where 

{F}  Is  a  vector  of  applied  loads 

<u>  is  a  vector  of  unknown  displacements  or,  as  in  the  present 
case,  displacement  Increments 

(K>  Is  the  global  stiffness  matrix,  in  the  present  case  it  is 
banded  and  positive  definite 
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INPUT 

FROM 

PERIFERAl/ 

DEVICE 


7 

!c? 


COMPUTATION  OCCURS  IN 
[1  PAIR  OF  BUFFERS  AS 
DATA  TRANSFER  OCCURS 
IN  OTHER  BUFFERS 


o 

O' 


OUTPUT  TO 
PERIFERAL 
DEVICE 


INPUT 


OUTPUT 


CORE  BUFFERS 


Initiate  read  of  data  into  (1) 

Test  for  completion  of  Step  a  data  transfer 
Initiate  read  of  data  into  (2) 

X,r:,3„ta(3)n  on  data  ,n  (,)  st°ri"9 

Initiate  write  of  data  in  (3) 

Test  for  completion  of  Step  c  data  transfer;  if  incomplete 
stop  computation  until  completed.  complete, 

Initiate  read  of  data  into  (1) 

Perform  computation  on  data  in  (2)  storing  results  in  (4) 
Test  for  completion  of  Step  e  data  transfer 
Initiate  write  of  data  in  (4) 

Test  for  completion  of  Step  g  data  transfer 
Loop  to  Step  c  until  all  data  has  been  processed 
Test  for  completion  of  final  buffer  operation 


FIGURE  5-22.  TYPICAL  OPERATION  OF  MULTI  BUFFERING  TECHNIQUE 
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Many  methods  of  obtaining  a  solution  to  Equation  5-5  .are  avai  lable.  A 
frequently  used  method  is  the  Choleski  Decomposition  Method.  Defining  1 

[K]  =  [L]  [D]  [L]T  (5-6) 


where 

[L]  is  a  lower  triangular  matrix  with  one.  on  the  diagonal  terms 
[D]  is  a  diagonal  matrix 


then  substituting  Equation  5-6  into  Equation  5-5 
{F}  -  IL]  [D]  [L]T  {u} 
and  defining 

{Z>  -  [D]  [L]T  {u} 

Equation  5"7  becomes 


(5-7) 


(5-8) 


m  -  id  {z>  ,5.,) 

There  are  many  algorithms  for  solving  Equations  5-8  and  5-9.  The 


th  i  s 

code 

are 

(Reference  5-4) . 

L.  . 
U 

= 

1 

D. 

J 

r  £  1 

KtJ  -S  L  im  J  ,m 

m=  l  j  j 

i  >  j  (5-10) 

D. 

J 

K.. 

JJ 

-  ^  D  L? 
mi  m  Jm 

(5-11) 

i  =  1 

Z. 

1 

Fi  “ 

£  Z- 
m=l 

(5-12) 

Z. 

n 

u. 

J 

s 

_1 

D. 

J 

-  5]  L  .  U 

m-J+1  mj  m 

(5-13) 
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where 

n  =  Total  number  of  equations 
i  ,j  =  Row  and  column  indices 


The  use  of  Equations  5-10  and  5“11  to  obtain  the  [L]  and  [D] 
matrices  in  the  most  general  case  of  a  block-by-block  solution  is  illustrated 
in  Figures  5"23  through  5“ 2 5 •  Figure  5 ” 23  shows  a  typical  banded  stiffness 
matrix  with  the  arrangement  of  the  terms  in  each  core  block.  Although  shown 
as  a  two-dimensional  array,  the  actual  locations  of  storage  may  be  consecutive. 
Figure  5-24  shows  a  method  of  assigning  core  storage  to  allow  the  double 
buffering  scheme  of  Section  5-2.1  to  be  used.  The  indicated  Scratch  Core 
area  may  be  of  any  size  (it  must  be  at  least  four  words  in  length)  and  is 
used  as  a  buffer  area  to  save  columns  of  the  reduced  stiffness  matrix  needed  for 
for  reduction,  i.e.,  decompos i t ion  of  fut  stiffness  blocks.  All  I/O  opera¬ 
tions  between  the  scratch  area  of  core  and  peripheral  devices  use  the  technique 
of  Section  5.2.1.  Thus  storage  and  retrieval  of  intermediate  data  required  for 
the  computation  of  [L]  and  [D]  may  in  the  optimum  case  be  performed  at 
core  speed.  Only  a  single  output  buffer  is  shown  in  Figure  5-24.  This  presents 
no  contradiction  to  the  double  buffering  scheme.  Figure  5-25  shows  that 
between  initiation  of  the  write  and  storage  of  new  data  in  the  output  buffer 
many  calculations  are  performed.  Tests  have  shown  that  this  calculation  time 
is  greater  than  the  data  transfer  time  thus  allowing  a  single  output  buffer 
to  be  used.  Figure  5“25  shows  the  sequence  of  operations  in  decomposing  the 
stiffness  matrix.  The  generation  of  {Z}  (from  Equation  5-12)  may  be  per¬ 
formed  at  the  same  time  [L]  and  [0]  are  generated. 


After  completing  decomposition  of  the  stiffness  matrix  and  generating 
the  {Z}  matrix,  the  displacements  {u}  are  computed  by  Equation  5-13- 
Figure  5~ 26  shows  the  sequence  of  operations  in  computing  the  {u}  matrix. 
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FIGURE  5-25.  SCHEMATIC  DIAGRAM  OF  STIFFNESS  MATRIX  DECOMPOSITION 
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FIGURE  5*26.  SCHEMATIC  DIAGRAM  OF  SOLUTION  VECTOR  EVALUATION 
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SECTION  6 

COMPARISON  WITH  CLOSED-FORM  ANALYTIC  SOLUTIONS 


6.1  SAMPLE  PROBLEMS 

To  Investigate  the  numerical  accuracy  of  the  present  computer  program 
and  to  determine  the  computer  time  required  to  solve  problems  of  various  size 
and  con?1exity,  several  example  problems  have  been  formulated  and  their 
numerical  solutions  cotyaredwlth  closed-form,  analytic  solutions.  One-,  two-, 
and  three-dimensional  linear  elastic  solutions  are  considered  as  well  as  one- 
and  two-dimensional  elastic/plastic  and  two-dimensional  visco-elastic  solu¬ 
tions.  These  are  listed  In  Table  6-'. 


Problem  1— Stress  Around  a  Circular  Hole 

The  geometry  of  Problem  No.  1  Is  illustrated  in  Figure  6-1.  The 
finite  element  mesh  Is  Illustrated  In  Figure  6-2.  The  solution  is  shown  in 
Figure  6-3  In  terms  of  principal  stresses  at  0  -  5-7  and  42.1°.  It  does  not 
depend  on  the  material  properties  of  the  plate. 


Problem  2— Stress  and  Displacement  In  an  Elastic 
""TfiTck-Wal  led  Cylinder  Under  Internal  Pressure 

The  geoiwtry  of  Problem  No.  2  is  Illustrated  In  Figure  6-4.  The 
finite  element  mesh  is  also  shown.  The  solution  is  shown  In  Figure  6-5  In 
term  of  radial  and  tangential  stresses. 


p r ob 1 em  3— Stress  In  an  Elastic,  Perfectly  Plastic 
Thlck-Ujl led  Cylinder  Under  Internal  Pressure 


The  geonvtry  of  Problem  No.  3  l»  the  same  as  that  of  Problem  No.  2 
and  Is  shown  In  Figure  6-4.  For  this  problem,  an  additional  material  property, 
the  Mises  yield  criterion,  Is  specified  as  follows: 


f  -  /jr-  a,  <  0 


(6-1) 


The  analytic  solution  is  shown  in  Figure  6-6. 
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TABLE  6-1.  PROBLEMS  SOLVED  BY  FINITE  ELEMENT 
AND  CLOSED-FORM  METHODS 


Nurtbc r 

Description 

Geomc t ry 

Material  Property 

Closed-Form 

Solution 

1 

Stress  around  a 
ci rcuiar  hole 

Two-dimensional 

(plane) 

Elastic 

Reference  6-1 

2 

Stress  in  a  thick 
wal led  cyi inder 
under  internal 
pressure 

One-d Imens iona 1 
(axi symmetric) 

Elastic 

Reference  6-2 

3 

Stress  in  a  thick 
walled  cyi inder 
under  internal 
pressure 

One-dimensional 
(axi symmetric) 

Elastic,  perfectly 
plastic 

Reference  6-2 

U 

Stress  in  a  rein¬ 
forced  thick 
wai led  cyi inder 

One-dimens ional 
(axi symmetric) 

Viscoelastic 

Reference  6-3 

5 

Stress  concentra¬ 
tion  around  a 
cylindrical  hole 
in  a  semi  inf ini te 
body 

Three- 

dimensional 

Elastic 

Reference  6-li 

•  31 
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SQUARE  PLATE  ( 20  IN.  X  20  IN.  I 


FIGURE  6-2.  FINITE  ELEMENT  MESH  FOR  STRESS  CONCENTRATION  AROUND 
CIRCULAR  HOLE 
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FIGURE  6-3. 


COMPARISON  BETWEEN  PRESENT  FINITE  ELEMENT  SOLUTION  *md  aiiaivt.i- 
SOLUnON  FOK  PROBLEM  I  STRESSES  AROUNOR  CIRUu!  Hole 


FIGURE  6-4.  PROBLEM  2  INFINITELY  LONG,  THICK  ELASTIC  CYLINDER  SUBJECT 
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Problem  ^--Stress  in  a  Visco-Elastic,  Reinforced 
Cylinder  Under  Internal  Pressure^ 

The  geometry  of  Problem  No.  4,  shown  in  Figure  6-7,  is  similar  to 
that  of  Problems  2  and  2-  The  main  difference  is  a  steel  reinforcing  ring 
around  the  outer  circumference.  The  material  of  the  cylinder  is  assumed  to  be 
governed  by  a  Maxwell -type  law  as  follows: 


°;j  ■  2G (e ! j)  exp  (- t/B) 


(6-2) 


where 

e!j  =  Component  of  deviatoric  strain  tensor 

G  =  Elastic  shear  modulus 

t  =  Time  in  units  of  B  =  n/G  where  n  =  viscoelastic  parameter) 

Volumetric  deformation  of  the  cylinder  and  all  deformations  of  the  reinforcing 
ring  are  assumed  to  be  linearly  elastic  inviscld.  The  variation  of  radial  and 
circumferential  stresses  are  shown  as  functions  of  radius  and  time  in  Figure  6-8. 

Problem  5 — Three-Dimensional  Stress  Concentration  Around  a 
Cylindrical  Hole  in  a  Semiinfinite  Elastic  Body 

The  geometry  of  Problem  No.  5  is  illustrated  in  Figure  6-9.  The 
stress  distribution  around  the  hole  near  the  stress  free  face  (Xj  -  X2  plane) 
is  appropriate  to  plane  stress,  while  in  the  interior  there  is  axial  stress 
along  the  axis  of  the  hole.  The  finite  element  mesh  is  shown  in  Figure  6-10. 

The  loading  condition  selected  for  this  example  is  uniaxial  stress  parallel 
to  the  X2-axis.  Thus  the  faces  parallel  to  the  X2~p  plane  are  stress  free  as 
is  one  face  parallel  to  the  X2*X^  plane.  The  finite  element  solution  is 
compared  with  the  analytic  solution  in  Figure  6-11. 
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FIGURE  6-6.  PROBLEM  4 — RE  I NFORCED  VISCOELASTIC  CYLINDER  SUBJECTED  TO  INTERNAL 
PRESSURE  (REFERENCE  6-4) 


PROBLEM  5--FINITE  ELEMENT  MESH  (ONLY  FIRST  THREE  LAYERS  ARE 
SHOWN  FOR  CLARITY.  COMPLETE  MESH  CONTAINS  EIGHT  LAYERS.) 
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(a)  AXIAL  DISPLACEMENT  ON  THE  STRESS  FREE  PLANE  p  ■  0  (NOTICE  THAT  THE 

roAu^T  SOLUTION  WHICH  ENCOMPASSES  A  FINITE  DOMAIN  DIVERGES 

ruc^cJm tc^d^VT ' C  S0LUTI0N’  WH,CH  CONSIDERS  AN  INFINITE  DOMAIN,  AS 
THE  FINITE  BOUNDARY  IS  APPROACHED) 


*  *  r/i 


(b)  AXIAL  STRESS  AT  EDGE  OF  THE  HOLE 

F,t5URE  PK0BLEM  5— THREE-DIMENSIONAL  STRESS  CONCENTRATION  (REFERENCE  6-4) 
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(d)  TANGENTIAL  STRESS  AT  EDGE  OF  HOLE  (SPRING  LINE) 


FIGURE  6-11.  (CONTINUED) 
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6.2  COMPUTING  TIME  REQUIRED  FOR  SOLUTION 

The  computing  time  regal  red  to  solve  Problem  5  end  several  small 
problems  on  a  Unlvac  1108  with  65,000  words  of  core  is  shown  In  Table  6-2. 
More  data  on  time  to  solve  various  sizes  of  problems  will  be  gathered  during 
the  second  phase  of  the  contract. 
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TABLE  6-2.  COMPUTING  TIME  REQUIRED  FOR  SOLUTIONS 
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APPENDIX  A 

LOGIC  OIAGRAMS 
FOR  MATERIAL  PROPERTIES 

Logic  Jlagrams  for  subroutines  In  the  material  property  package 
are  shown  ielow.  Subroutine  CONECT  connects  the  package  to  the  main  program 
Subroutine  ELPL  controls  all  the  othe-  subroutines. 
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FIGURE  A-2.  SUBROUTINE  ELPL  CONTROLS  MATERIAL  PROPERTY  SUBROUTINES  TRANSFER' 
AXES  FOR  ANISOTROPIC  MATERIALS,  PERFORMS  TESTS  FOR  STRAIN 
SPLITTING  AND  ITERATION. 
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FIGURE  A-l».  SUBROUTINE  PLAST— THIS  SUBROUTINE  MODIFIES  C  MATRIX 
GENERATED  BY  SUBROUTINE  ELAST  IN  ORDER  TO  ACCOUNT  FOR 
PLASTICITY.  THE  BASIC  QUANTITIES  FOR  THIS  MODIFICATION 
ARE  DERIVATIVES  OF  THE  YIELD  FUNCTION  WITH  RESPECT  TO 
THEIR  ARGUMENTS  AND  ARE  COMPUTED  IN  YLDFUN. 
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SUBROUTINE  TSTYLD — COMPARES  CURRENT  STATE  OF  STRF<;<;  uitu  vicim 
CRITERIA  TO  DETERMINE  WHETHER  YIELDING  OCCURS  Y  ELD 
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START 


FIGURE  A-9.  SUBROUTINE  VPLAST— THIS  SUBROUTINE  IS  USED  FOR  A  VISCOPLASTIC 
MATERIAL.  IT  CALCULATES  THE  STATIC  YIELD  CRITERION,  DETERMINES 
WHETHER  YIELDING  OCCURS  AND,  IF  YIELDING  OCCURS,  MODIFIES  THE 
STRESSES  TO  ACCOUNT  FOR  VISCOPLASTIC  STRESS  RELAXATION. 
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FIGURE  A-10.  SUBROUTINE  COMPUT— COMPUT  CALCULATES  STRESS  INCREMENTS  AND  ADDS 

THEM  TO  OLD  STRESSES  TO  OBTAIN  NEW  STRESSES.  STRAINS  CORRESPONDING 
TO  ZERO  STRESSES  (PLANE  STRESS,  UNIAXIAL  STRESS) ARE  COMPUTED. 
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